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ABSTRACT 



New concepts in aerospace travel have renewed interest 
in modeling and compensating for the effects of upper 
atmospheric drag. In particular, the SDI constellation 
requires strict orbital element maintenance. This thesis is 
a qualitative verification of a propellant longevity model 
for low-altitude earth orbit satellites doing intrack micro- 
thrusting to overcome atmospheric drag. The original model 
was developed at the Naval Surface Weapons Center. 
Pertinent orbital mechanics and atmospheric concepts are 
reviewed. The model and its computer program are described. 
The results of trend and sensitivity analysis reasonableness 
tests are presented. Finally suggestions for use are made. 
Many plots of mission life predictions are presented. The 
model computer program and sample input and output are also 
included. 
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I. 



INTRODUCTION 



All current low-earth-orbit satellites are simply placed 
in orbit and then tracked as their orbits decay deeper into 
the atmosphere. This orbital decay is due to forces on the 
satellite produced by atmospheric drag and the non- 
spherical earth. 

Modified exponential-t^e models for the atmosphere 
produce sufficiently accurate drag predictions to track the 
satellites. 

New concepts in scientific and military uses of space 
require greater accuracy in predicting atmospheric drag 
effects then produced by the modified exponential models. 
These new uses include the Aero-assisted Orbital Transfer 
Vehicle (AOTV) , orbiting space stations, the Trans 
Atmospheric Vehicle (TAV) , and the strategic defense 
initiative satellite constellation. 

The environment of all of these vehicles is called the 
thermosphere, an atmospheric layer studied in depth only 
recently. An accurate empirical model of the thermosphere 
has also been produced. 

One idea under investigation at the Naval Surface 
Weapons Center is the placement of many low-earth-orbital 
defensive satellites in precisely controlled and maintained 
orbits. This set of defensive satellites is termed a 
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constellation. All low earth-orbital-satellites drift from 
their initial orbits due to perturbating accelerations 
caused primarily by the aspherical earth and atmospheric 
drag. In order to be functional, the orbits of the 
satellites in an SDI constellation must be exactly 
maintained. The drift of a satellite must be cancelled out 
by using an onboard propulsion system to produce forces to 
offset the perturbating accelerations. 

One model, proposed by Dr. Dan Parks of the Naval 
Surface Weapons Center, calculates the mission lifetime 
expected of a satellite doing intrack micro-thrusting to 
overcome atmospheric drag effects. The end of mission life 
corresponds to depletion of all the maneuvering fuel. This 
model is here-after referred to as the propellant longevity 
model. The propellant longevity model incorporates a Bessel 
function expansion for the fuel-mass decrement on each 
revolution. However, the model does not calculate the fuel 
necessary to cancel out orbital drift due to the aspherical 
earth . 

This thesis gives a verification of the low-altitude 
earth satellite propellant longevity prediction model 
developed by Dr. Parks. The theory for the model 
development itself is presented in Naval Surface Weapons 
Center Technical Report (NSWC TR) 83-243 (Parks [Ref. 1]). 
Computer code for the model was developed at NSWC under the 
direction of Dr. Parks. In addition to the computer code of 
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the model, several other NSWC publications with applications 
to the model were forwarded to the Naval Postgraduate 
School . 

The Thesis is organized as follows. Chapter I presents 
a general introduction to the propellant longevity model and 
discusses the changes made to the computer code for use at 
the Naval Postgraduate School. A basic description of the 
orbital mechanics necessary to operate the model is provided 
in Chapter II. Chapter II also presents an overview of the 
atmospheric physics necessary to understand the density and 
compositional variations in the thermospheric layer of the 
thermosphere. A description of the propellant longevity 
model is given in Chapter III. Instructions for operation of 
the computer code are also included there. Chapter IV 
contains the results of the verification process sensitivity 
analysis and reasonableness tests. A summary is given in 
Chapter V. Appendix A contains the propellant longevity 
Fortran computer code for the IBM 3033. Appendix B shows 
input sets and their corresponding outputs. Graphical 
presentation of the tabular results is included. 

The coding of the model by NSWC was completed in Fortran 
for the Control Data Corporation (CDC) 6700 computer. 
Because of software and hardware differences, extensive 
modifications to the model coding were required to install 
the program as an operational model on the IBM 3033 computer 
at the Naval Postgraduate School. To optimize the utility 
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of the model, some adjustments were also made to the 
original computer coding of the propellant longevity model. 
Great care was taken when modifying the code to preserve the 
model theory as contained in NSWC TR 83-243. 

The alterations to the computer program fall into two 
categories: changes to transport and implement the model on 
our IBM 3033, and changes made to mainstream the. model 
verification process. Changes to transport the model from 
the CDC environment to the IBM environment involved some 
logical revision and syntax, hardware word-size 
cons iderat ions . 

The hardware word-size problems were extensive. The CDC 
6700 uses a 60-bit word size for single precision 
calculations, whereas IBM's standard is a 30-bit word. The 
propellant longevity model requires 60 bit precision. Thus 
all real numbers, library functions, and input constants 
were promoted to double precision accuracy. Another 
precision-related problem involved the use of fractional 
exponentiation. All one-half exponents were changed to call 
the SQRT function to increase precision and speed. 

Syntactical changes were mostly related to the software 
environment transportation from CDC Fortran to IBM Fortran 
77. Library functions calls, and input and output format 
statement changes, were the bulk of the syntactical changes 
made. (Examples of library function changes are the ARQTAN 
function changed to the ATAN2 function, the BESI Bessel 
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function call changed to the MMBSIR Bessel function call and 
so forth) . 

Examples of changes in the logic to in order to 
transport and implement the model code include using the 
Fortran 77 block IF-THEN-ELSE statement for the CDC Fortran 
IF statement and the use of DIMENSION statements in each 
subroutine that is processing a particular array. 

The second category of changes to the program addressed 
mainstreaming the model verification process. The model 
code as down-loaded from the NSWC tape contained little 
code-internal docximentation . Looking towards future 
implementation and maintenance of the model code, extensive 
documentation was added as code internal comments. Most of 
this documentation was based on the research literature of 
NSWC. Extensive reference to the current publications is 
made in the code documentation itself. 

Two changes were made in the way the code solves the 
model in order to implement the theory with greater 
precision. The first change turned on the subroutine THRST 
which was designed to allow for preplanned changes to the 
semi-major axis of the orbit at the expense of onboard, 
maneuvering fuel. (The subroutine code down-loaded from 
tape did not return the fuel consumed for the orbit adjust, 
record the resultant change in the semi-major axis, or 
receive the correct intended adjust parameter. All these 
functions are performed in the IBM 3033 version.) 
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A second change in the code alters the way in which the 
satellite velocity relative to the atmosphere is computed. 
In the original code this parameter had to be entered as 
part of the input file and was assumed to be constant. 
Since this relative velocity is actually a function of other 
transient orbital parameters, recomputation is now 
automatically accomplished for each satellite revolution 
based on current values of the changing orbital parameters. 
In place of the satellite-to-atmospheric-relative-velocity 
input, the user now inputs the atmospheric-to-earth relative 
rotation rates. The physics associated with this change is 
discussed later on, in the section describing the 
atmospheric model. 

A few other changes to the propellant longevity code 
customize the output, or state the input element initial 
conditions. Several control inputs now allow for a range of 
model predictions. These are explained in Chapter III. 

The actual model verification is accomplished in three 
steps. First does the model address the original problem? 
Second, does the model meet the. criteria of common sense? 
Finally, how well does the model solution compare when 
tested against real-world data? 

No real-world data sets exist for testing the propellant 
longevity model. The only low-earth-orbit satellite having 
a propulsion system to maintain orbital parameters is highly 
classified under the auspices of the United States Air 
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Force. The data from this program is not available in an 
unclassified state. In lieu of actual data, several 
additions to the program code were made to assist in the 
model verification process. The major adjustment is the 
addition of a feature to allow for a sensitivity analysis. 
The subroutine SNALYS, allows the user to linearly scan any 
combination of the 20 input elements, from a user specified 
minimum to any maximum, with user specified granularity. 
The output of the sensitivity analysis can be descriptive 
and tabular, or graphical in nature. The SNALYS subroutine 
allows for a plot of model prediction against the entire 
range of the input elements. Results of input element 
sensitivity analysis are presented in the Chapter IV. The 
sensitivity analysis feature allows for optimum orbital- 
element and satellite-ballistic-coefficient design. 

Thought was given to possible model installation on 
other computers. The model is large and requires about 20 
minutes of dedicated central processor time to complete an 
average run. As such the model is rather limited to 
mainframe computers. Initial calculations based on 
processor time indicate an IBM AT microcomputer with the 
requisite 80287 co-processor and additional RAM memory might 
require about ten hours to complete an average run. 

This chapter presented an overview and motivation for 
the propellant longevity .model. Changes made to the 
computer program at the Naval Postgraduate School were also 
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discussed. The next chapter presents the necessary 
background in orbital mechanics and atmospheric concepts to 
operate and understand the model. 
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II 



BACKGROUND 



This chapter presents the orbital mechanics and 
atmospheric concepts required for the understanding and 
operation of the propellant longevity model. The orbital 
mechanics presented in the first section discusses only the 
orbital element sets, concepts, and coordinate systems 
specific to the model. In the second section of this 
chapter. Thermosphere specific atmospheric concepts and 
cycles are presented. 

A. ORBITAL MECHANICS 

This section covers the orbital mechanics for operating 
the mass decrement model. If the coordinate and orbital 
element systems discussed are familiar, this basic orbital 
mechanics section can be skipped. Orbital mechanics 
references are cited in the text. 

Early celestial mechanics models were based on 
increasingly accurate observations. As observations and 
record keeping improved so did the empirical model fit to 
the observations. The model that fathered modern orbital 
mechanics was the empirical model derived by Johannes 
Kepler, 1571-1630. In 1619 Kepler formulated three laws 
(models) of planetary motion. These models were based on 
life long observations recorded by Tycho Brahe, 1564-1601. 
Kepler's models were formulated as follows; 
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1. The orbit of each planet is an ellipse with the sun at 
one focus. 

2. The line joining a planet to the sun sweeps out equal 
areas in equal times. 

3 . The square of the period of a planet is proportional 
to the cube of its mean distance from the sun. 

Later, Sir Isaac Newton, 1643-1727 validated Kepler's 
laws by establishing their equivalence with his law of 
universal gravitational attraction. 

Newton's model for gravitational attraction states that 
two bodies attract with a force that is proportional to the 
product of their masses and inversely proportional to the 
square of the distance between them. Symbolically, the 
force on the mass m due to the mass M is given by: 



Fg 



GMm r 
2 r 



r 



( 1 ) 



where G is the universal gravitational constant and r is the 
vector from the center of mass of M to the center of mass m. 
r is the magnitude of the distance vector r. 

A body experiences an acceleration due to the resultant 
of all the forces acting on the body. If only a 
gravitational force acts on the body, then the acceleration 
of the body is given by: 



r" = a 



G(M+m) - 

3 ^ 

r 



( 2 ) 
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In this model, r" is the acceleration vector on m caused by 
gravitational attraction of M on m. The remaining variables 
are defined as in the previous model. If M is much greater 
then m (as in an earth satellite system) , G(M+m) is very 
nearly GM. Another constant can then be defined called the 
gravitational potential parameter u: 

u = Gm (3) 

The gravitational acceleration of a satellite in earth 
orbit can then be expressed as: 

r" + u/r^ * r = 0 (4) 

The gravitational field is conservative. This means 
that an energy constant of motion can be defined for a body 
moving in a gravitational field. For a two body, an- 
atmospheric, elliptic orbit system the specific mechanical 
energy is: 

^ = vt^/2 - u/r (5) 

In this model ^ is the specific mechanical energy, v-j- is 
the satellite tangential speed, and u and r are as 
previously defined. Note that E is the sum of the kinetic 
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energy per unit mass of the satellite and its potential 
energy per unit mass. 

The model, (5), requires satellites at specific radii to 
have specific tangential speeds. As r increases, v must 
decrease to maintain C constant. This orbital energy model 
is important when modelling the effects of atmospheric drag 
which is the major topic of the next section. 

Newton's inverse square law of gravitational force 
requires gravitational field orbits to be elliptical in the 
shape. The mathematics of the ellipse is important in the 
description of orbits and Figure 1 presents some of the 
important parameters. 

The parameters shown in Figure 1 are defined as follows; 
a: the semi-major axis, 

b: the semi-minor axis, 

e: the eccentricity where 

e = SQRT [1 - (b/a)2]. 



C: 



0 : 



p: 



r : 



v: 



rp: 



ra: 



the ellipse center. 

the ellipse focus (earth center) . 

the semi-latus rectum. 

the distance from a position on the ellipse 
the focus 0. 

the true anomaly. 

the radius at perigee (in an earth system) . 
the radius at apogee (in an earth system). 
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Figure 1. The Elliptical Orbit 

A: the apogee (point farthest from the focus) . 

P: the perigee (point closest to the focus. 

No real orbits occur in just two dimensions. To describe 
real orbits a third dimension must be added. A suitable 
reference must also be chosen. One of the most convenient 
reference frames is the geocentric reference depicted in 
Figure 2. In this system a fixed nonrotating earth forms 
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the center of the celestial sphere. The earth's projected 
north pole, south pole, and equator become the north 
celestial pole (NCP) , the south celestial pole (SCP) , and 
celestial equator respectively. 




Figure 2 . The Geocentric Celestial Sphere 

To account for the apparent motion of the celestial 
sphere as viewed from the earth, a reference point on the 
celestial sphere must be chosen. Traditionally this has 
been the First Point of Aries, describing the constellation 
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where the vernal equinox occurred when the celestial 
reference was defined. The vernal equinox does not now 
occur in the constellation Aries. Tradition still uses the 
first point of Aries as the Celestial Meridian (CM) of 
reference. Figure 3 illustrates the vernal equinox 
reference. The vernal equinox is defined as the time when 
the sun is directly over the equator at noon in the spring. 




Figure 3 . The Vernal Equinox 
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Figure 4 is a representation of an orbit in real space 
using the geocentric system. The points where the 
projection of the satellites orbit on the celestial sphere 
crosses the celestial equator are called the nodes of the 
orbit. The ascending node (N) is where the projected orbit 
crosses the celestial equator heading for the north 
celestial pole. The descending node is where the projected 
orbit crosses the celestial equator heading for the south 
celestial pole. The angular measurement from the reference 
celestial meridian to the ascending node along the celestial 
equator is called the longitude of the ascending node, and 
usually given the symbol . The angular measurement from 
the ascending node to the point of perigee (P) along the 
projected orbit is called the argument of perigee, 
symbolized by co . Both the longitude of the ascending node 
and the argument of perigee range from 0 to 360 degrees. 
The angle from the plane of the celestial equator to the 
plane of the projected orbit is called the inclination and 
denoted by i. Values for the inclination range between 0 
and 180 degrees. 

Inclinations of 0 and 180 degrees are termed equatorial. 
In equatorial orbits the longitude of the ascending node is 
not defined. Orbits with an inclination of 90 degrees are 
called polar. Prograde orbits are those with an inclination 
of between 0 and 90 degrees and retrograde orbits have an 
inclination range of 90 to 180 degrees. 
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NCP 




Figure 4 . A Geocentric Orbital System 

Historically, it was difficult to tell much about the 
orbit specific parameters from the earth vantage point. 
However, a value that could be measured was the time for one 
complete orbit. If the radius vector moved uniformly along 
the projected orbit with the average rate of mean motion n 
(n is also referred to as the mean angular velocity) at an 
epoch T, it would be at a position M called the mean 
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anomaly. The mean anomaly has a range of 0 to 360 degrees. 
It is measured from the perigee along the projected orbit. 
The mean anomaly corresponds to positions on the orbital 
ellipse. An eccentric anomaly (E ) , is defined as the angle 
between the projected position and the point of perigee 
subtending at the center of the ellipse. Figure 5 
illustrates eccentric anomaly. 




Figure 5. Orbital Anomalies 
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The difference between the eccentric anomaly and the mean 
anomaly is important. The mean anomaly is a mathematical 
concept corresponding to an imaginary position of a 
satellite past the perigee if it moved at an average angular 
velocity. The eccentric anomaly is the actual projected 
angular position of a real satellite on the celestial 
sphere. The eccentric and mean anomalies are equal in 
circular orbits (e = 0) . The eccentric and mean anomalies 
are related by Kepler's equation: 

M = E - e*sin E (6) 

The solution to model (6) has been the source of much 
investigation. By a solution it is ' meant for a given M, 
determine e and E. In 1817 Friedrich Wilhelm Bessel, 1784- 
1846 invented Bessel functions to provide an infinite series 
approximation for the solution of the model (6) . Newton 
developed a converging algorithm for its solution. The mass 
decrement submodel of the propellant longevity model uses 
Newton's method for the solution of Kepler's equation. 

All anomalies are concerned with the measurement of the 
position of a body in orbit. All anomalies are measured 
angularly from the point of perigee in the orbital plane, 
and vanish at perigee. 

When e and E have been found from Kepler's equation (6), 
then using the equations for an ellipse in Taff [Ref. 2:p. 
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67], the remaining elliptical parameters can be solved for 
when any one of a, b, rp, or ra is known. 

The minimum number of parameters needed to describe an 
orbit is called an orbital element set. Five parameters are 
needed to describe an orbit. One more is needed to position 
an object in the orbit. Many orbital element sets have been 
used and are in use. The classical orbital elements are: 
semi-major axis (a) 
eccentricity (e) 
inclination (i) 

longitude of the ascending node (f2) 

argument of perigee (oj) 

time of perigee passage or epoch (T) . 

The Keplerian set replaces the time of perigee passage with 
the mean anomaly (M) . The orbital element set used in the 
mass decrement submodel is the Kozai mean element set, very 
similar to the Keplerian set. The name "Kozai mean element 
set" refers to the Kozai perturbation approximations used to 
solve for the Keplerian orbital element set. The 

perturbation theory approximation of Y. Kozai is discussed 
later. 

Up to now all orbital considerations have been limited 
to the classical two-body problem. The two body problem 
limits the forces on the orbiting body to the gravitational 
force. The two-body problem also requires all masses be 
concentrated at a single point. A real earth satellite 
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experiences many forces resulting in many discrete 
accelerations. The theory of perturbations was developed to 
account for the motion of real satellites. 

In perturbation theory the original orbit, with only the 
point -mass gravitational acceleration, is called the 
osculating orbit. Orbital elements describing the 
osculating orbit are the osculating element set. The 
osculating orbit is then adjusted periodically to account 
for the action of all the nonpoint-mass gravitational 
accelerations . 

Accelerations result from forces. Examples of 
perturbating forces on an earth satellite are: oblate earth 
gravitational, atmospheric drag, sun gravitational, moon 
gravitational, solar radiation pressure, and satellite motor 
thrust. 

To remain in a closed orbit, the point -mass 
gravitational force must be the dominant force on the 
satellite. For low earth orbit satellites, the oblate earth 
gravitational force is the next greatest force and is termed 
the dominant pertiirbation. Satellites below 600 kilometers 
in altitude experience an atmospheric drag force of the same 
order of magnitude as the oblate earth force. The 
atmospheric drag and oblate earth effects are so large for 
low earth orbit satellites that the effects of the other 
perturbations (except motor thrust) are . lost in the 
background error. 
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Central to the problem of perturbation theory is the 
that the three body problem has never been solved 
mathematically in general closed form. Most modern 
approximations involve power series or Legendre polynomial 
expansions. The state of the art in orbit perturbation 
prediction only uses only the first-order terms of these 
expansions, termed first-order effects. An exception is a 
theory advanced in 1959 by Y. Kozai. His approximations are 
computationally efficient and allow for the use of more of 
the oblate earth effects (called J2, J3, and J4 terms and 
discussed later). Kozai 's theory is flawed in the higher- 
order terms due to the introduction of mathematical 
singularities. Nevertheless, proponents for its use argue 
that it does work and predicts the oblate earth effects to a 
reasonable degree of accuracy. The mass decrement submodel 
of the propellant longevity model presented in this thesis 
uses Kozai 's theory. 

The gravitational potential parameter u was presented as 
singularly determining the gravitational potential in the 
earlier development of the gravitational acceleration model. 
In earth systems, the gravitational potential is a complex 
function of the earth's mass distribution. 

The earth is not a point mass. Neither is it a 
homogenous sphere. It's true mass distribution is complex. 
The real distribution of matter on the earth forms a body 
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termed the geoid. This geoid determines exactly the 
gravitational potential function. 

The determination of the geoid is based on ground and 
satellite measurements. Many military programs depend on 
the precise position of a point on the surface of the earth. 
The military has led the way in geoid determination in the 
past two decades. 

Geoidal studies reveal important facts about the 
aspherical shape of the earth. A fluid spinning sphere 
should exhibit a bulge at the equator and flattening at the 
poles. This is termed oblateness. The oblate deviation is 
the dominant deviation for the earth as a sphere. 
Measurements of the earth’s oblateness indicate that it is 
more oblate than the current spin rate could produce by 
itself. This excess oblateness indicates two factors; the 
earth's core is more solid than previously thought, and the 
earth's spin has slowed. 

The earth deviates in other ways from being a sphere. 
Since all deviations can be expressed in terms of sines and 
cosines they are referred to as spherical harmonics. 
Harmonics are divided into categories. Zonal harmonics are 
symmetric about the earth ' s spin axis and are not longitude 
dependent. Sectorial harmonics are only longitude 
dependent. Harmonics dependent on both longitude and 
latitude are termed tesseral. Odd numbered harmonics are 
antisymmetric about the equatorial plane; even numbered 
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harmonics are symmetric about the equatorial plane. 
Harmonics are also classified by number. The number refers 
to the subscript of the J coefficient (defined later) . 

Bate et al., [Ref. 3], offers a simplified development 
of the gravitational potential function. If the 
gravitational potential function $ of a geoid is known, the 
acceleration due to this gravitational function is modeled 
by: 



r" = V$ (7) 

One common expression for the gravitational potential or 
geoid function $ is: 

00 ‘ 

$ = -II - y J (— )^ P sin L] . (8) 

r n r n 

In this model u is the gravitational potential parameter, Jn 
is an experimentally determined J coefficient, r^ is the 
radius of the earth, Pn is a Legendre polynomial of order n, 
L is the geocentric latitude, and r is the radius of the 
satellites position. Taking the gradient of the potential 
function yields an expression for the acceleration due to 
the gravitational potential. Bate et al., [Ref. 3:pp. 419- 
423] details the resultant potential induced acceleration to 
the first seven terms. 
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Recent studies have determined the geoid function out 
through the J44 term. Values used by the mass decrement 
model by the subroutine GEOP are 1082.76*10“^, -2.55*10“®, 
and 1.56*10“®, for the J2, J3, and J4 terms respectively. 
The J2 and J4 terms are even numbered, and refer to even 
harmonics symmetric about the equator. The J3 term is odd 
numbered, and refers to odd harmonics which are 
antisymmetric about the equator. The J2, J3, and J4 terms 
are all zonal. 

Current perturbation theory expresses $ • to the J2 
term. The theory advanced by Kozai in 1959 appears to use 
the J2, J3, and J4 coefficients in a first-order 
perturbation theory approximation (because singularities 
occur in terms higher-order than J2) . Taff [Ref. 2:pp. 332- 
342] takes strong exception to the theory. 

The subroutine GEOP in the propellant longevity model 
uses Kozai 's theory to calculate the geoidal effects on 
orbital parameters through the first three zonal 
coefficients (J2, J3, and J4 terms). Their effect is 
calculated once each orbital revolution, and the Kozai mean 
element set is adjusted accordingly. While using the model, 
it is important to understand that as the altitude of the 
satellite decreases. The aspherical earth perturbations 
become much greater. The orbital element drift due to the 
aspherical earth is thus much faster in lower altitude 
orbits. 
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In siommary, the computer program of the propellant 
longevity model uses three orbital element sets in a 
geocentric reference. The input orbital element set is the 
Brouwer mean element set of the Naval Space Surveillance 
System (NAVSPASUR) ^ system. Chapter III and Appendix A 
explain the input elements in detail in the input 
description. The Brouwer mean elements are first converted 
into a Brouwer osculating set. These osculating orbital 
elements are then transformed into a Keplerian-like Kozai 
mean set. Half of the computer program of the propellant 
longevity model transforms the Brouwer mean element set into 
the Kozai mean set. These element set transformations are 
used to prepare real-world orbital element sets for use in 
the propellant longevity predictions. Once the computer 
program has solved for the Kozai mean element set all 
perturbation effects result in adjustments to the Kozai mean 
set. Adjustments are made once per revolution. 

B. THE ATMOSPHERE MODEL 

Central to the validation of any model is an 
understanding of the system to be modeled. The model under 
investigation determines the propellant longevity for low 
earth orbital satellites doing in track micro-thrusting to 

^The NAVSPASUR system publishes a semiannual report on 
the 1000+ orbiting objects they track in The Satellite 
Situation Summary . The orbits are described in a five card, 
80 colximn, Fortran style input format. The report is 
available on request from the Commanding Officer, Naval 
Space Surveillance System, Dahlgren, Virginia. 
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overcome the effects of atmospheric drag. The fuel required 
by a particular rocket motor on a satellite to perform a 
known amount of work can be computed from the rocket 
equation. The problem occurs in computing the amount of 
work required. 

To compute the required work by the satellite propulsion 
system to overcome atmospheric drag requires knowledge of 
the atmospheric environment of the satellite. For low earth 
satellites the atmospheric environment is the thermosphere. 
The physics of the thermosphere is discussed in this 
section. Before discussing the thermosphere, some 
preliminary thermospheric properties are presented. 

For low earth orbital satellites, the drag- induced 
change in the orbital elements is significant. The oblate 
earth (referred to as the "J2 term”) orbital perturbation is 
normally the dominant perturbation term. The perturbations 
due to drag are of the same order of magnitude as those due 
to the oblate earth for orbits below 600 kilometers. 

The drag force caused by the satellite moving through 
the atmosphere always acts against the motion of the 
satellite. Neglecting the lift effects due to the 
satellite's shape and attitude, the drag force causes work 
to be done on the atmosphere by the satellite. 

From Bate e al., [Ref. 3:pp. 15-16] a satellite in orbit 
has a constant specific mechanical energy C given by 
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where C is the mechanical energy per unit mass, v-j- the 
satellite speed, u the gravitational potential constant, and 
r the radius of the orbit. 

The above model says that the mechanical energy per unit 
mass of a body in an an-atmospheric , two body, circular 
orbit is the sum of its kinetic energy per unit mass and 
gravitational potential energy per unit mass. 

Conservation of energy requires that the work done on 
the atmosphere by the satellite must reduce the orbital 
energy of the satellite. The differential element of work 
done on the atmosphere by the satellite results in a 
differential reduction in the specific kinetic energy of the 
satellite. Assuming the specific mechanical energy is a 
constant, the radius, r, of the orbit must therefore 
decrease. As the radius decreases, the velocity tends to 
increase in order to maintain the an-atmospheric 
conservation of specific mechanical energy. 

In non-circular (real elliptical) orbits, the 
atmospheric effects are greatest at perigee (Figure 1) . 
Orbital mechanics requires a differential reduction in 
orbital velocity at perigee to result in a reduction in the 
semi-major axis and orbital eccentricity. The major effect 
of transatmospheric satellite motion is to reduce the semi- 
major axis and the eccentricity of the elliptical orbit. 
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Taff [Ref. 2:pp. 352-353] develops the mathematical 
equations for the drag-induced changes in the semi-major 
axis, eccentricity, argument of perigee (Figure 4) , and mean 
anomaly (Figure 5) . The changes in the argument of perigee 
and the mean anomaly are shown to be periodic. The net 
effect on the orbit due to atmospheric drag is to 
circularize the orbit closer to the earth, with little 
change to the other orbital elements. 

An accepted model for the magnitude of the drag force as 
presented in Taff [Ref. 2;p. 352] is: 



Fd 



A pC^ 

2m ^rel 



( 10 ) 



In this equation, F^ is the drag force magnitude, A the 
effective satellite cross-sectional area, Cq the drag 
coefficient, m the satellite mass, p the atmospheric 
density, and v^g]^ the velocity of the satellite relative to 
the atmosphere. 

The upper atmosphere is dynamic. The changes that must 
be modeled to achieve a valid prediction for the drag force 
are those that effect the drag force factors. The factors 
in the drag force model that are dependent on the atmosphere 
are density, a drag coefficient (dependent on satellite 
physical 'and fluid-flow characteristics) , and the relative 
velocity of the satellite in the atmosphere. 
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To model the atmospheric effects on density, drag 
coefficient, and relative velocity, some background in upper 
atmosphere physics is needed. The descriptive presentation 
here is carefully developed to be relative to the above 
three drag force dependent factors. A more comprehensive 
treatment of atmospheric physics can be found in a work of 
Fleagle [Ref. 4]. An excellent technical treatment of 
thermospheric dynamics can be found in the works of Jacchia 
[Refs. 5,6,7]. 

Jacchia formulated the most widely used thermospheric 
model currently employed by professional orbital analysts 
(known as the Jacchia J60 model or simply J60) . His latest 
comprehensive empirical thermospheric model (Jacchia J77 or 
J77) is many orbital analysts' preferred refinement for 
modeling the atmosphere. 

The Jacchia 311 model is highly tabular, uses 
considerable computer time, and requires an extensive input 
file. However, these negative aspects are balanced by the 
orders-of-magnitude increase in . accuracy of density 
predictions. Jacchia [Ref. 7;p. 2] presents a table of 
residuals for predicted versus observed densities for 
various satellites, as well as time intervals for the 
Jacchia J77 atmosphere model. All density residuals, are 
much smaller than those produced by exponential models in 
the same time intervals. When extremely accurate density 
predictions are required in some submodel, the use of the 
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Jacchia J77 atmosphere model may well be worth the 
additional computer time and expense. 

A simple model of the atmosphere considers it as a 
homogeneous mass of constant density. A refinement is made 
by modeling the density as decreasing exponentially with 
altitude. The exponential model predicts pressure and 
density will fall one tenth for every ten kilometer increase 
in altitude. Prior to 1960 this exponential model was used 
extensively for upper air. calculations involving meteor 
trails. With the advent of low earth orbital satellites, a 
greater degree of accuracy in predicting atmospheric effects 
is required. The predicted satellite decays from the 
exponential model are very different from the observations. 
The Jacchia J60 atmosphere model was formulated to improve 
decay predictions. It was based on measurements form high 
altitude balloons, sounding rockets, and ground stations. 

The atmosphere is divided roughly into concentric 
shells, each with specific properties. The properties used 
in this subdivision are after Fleagle. [Ref. 4] 

From the early days of atmospheric studies, temperature 
and moisture have been recorded. The division system 
described by Fleagle [Ref. 4] is based primarily on a height 
versus temperature profile. 

The height versus temperature profile is graphed in 
Figure 6. The temperatures and altitudes are those of the 
1975 U.S. standard atmosphere. This temperature behavior 
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Figure 6. Atmospheric Temperature Profile 

approximates general conditions (large transient variations 
have been observed) . Figure 6 also shows the boundaries 
between concentric shells, called pauses. The following 
description of the atmospheric layers is presented in 
Fleagle [Ref. 4;pp: 79-84]. 

The region of the atmosphere closest to the earth is the 
troposphere and contains eighty percent of the atmospheric 
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mass. This region is characterized by near linearly 
decreasing temperature with altitude. The dynamics are due 
to earth surface energy transitions. All of the earth's 
weather occurs in this region. The height of the 
troposphere is variable and depends on both latitude and 
season.’ Along the equator the maximum altitude is 16 to 18 
kilometers. At the poles the troposphere is highly 
seasonal, being nearly absent during the winter months and 
extending to eight to ten kilometers in the summer. Capping 
the troposphere is a region of constant temperature with 
altitude, called the tropopause . 

Next in the series of atmospheric layers is the 
stratosphere . The stratosphere is characterized by linearly 
increasing temperature with height. The physical dynamics 
are due mainly to the absorption of solar ultra-violet 
radiation by the ozone. The maximum altitude of the 
stratosphere is about 50 kilometers where it is capped by an 
area of local maximum temperature termed the stratopause . 

Above the stratopause is the mesosphere . This 
atmospheric layer exhibits linearly decreasing temperature 
with height. The energy physics of the mesosphere are not 
well understood. Bulk mixing of the gasses still dominates 
diffusion as thermodynamic processes. Any molecular 
dissociation of atmospheric molecules produced by radiation 
absorption is rapidly consumed by chemical recombination of 
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these dissociated species. The mesosphere is topped by the 
turbopause at an altitude of about 90 kilometers. 

Above the turbopause the gross physics of the atmosphere 
changes drastically. This region is the environment of all 
current low earth orbit satellites. Termed the 
thermosphere . it extends from the turbopause to an altitude 
of 1000 to 2500 kilometers. Here bulk molecular mixing is 
dominated by diffusion. Lighter gasses diffuse upwards and 
heavier ones remain at the lower altitudes. Disassociation 
of molecular species becomes increasingly important, this 
results in an increase of individual atomic species. 
Standard flow equations (Navier-Stokes) are extremely 
difficult to apply, and physical descriptions are more 
discrete (quantum mechanical) . Flight by flow-supplied lift 
is not possible. The standard methods of computing drag 
coefficients do not apply. The ideal gas approximation is 
also not valid (P*V = n*R*T does not apply) . Above the 
turbopause, density is an increasing function of 
temperature. Only about one millionth of the atmosphere 
remains above the tropopause. 

The stratification of gaseous components in the 
thermosphere is significant. The drag coefficient has been 
shown to be composition dependent. Thus, as the gaseous 
composition changes the drag coefficient must be adjusted to 
reflect the observed drag-induced orbital decays (Jacchia 
[Ref. 7]). Figure 7 shows the dominant atmospheric species 
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Figure 7. Atmospheric Composition 

by altitude with the accepted drag coefficients. The drag 
coefficients are predictions obtained analyzing the drag 
coefficient needed to produce the observed drag effects on 
orbital decay for a satellite at a particular altitude. 
Thus the drag coefficient is a proportionality constant in 
the drag force model. 
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For altitudes below 110 kilometers, molecular nitrogen 
is the main atmospheric gas. Above this region the 
dissociation of the molecular oxygen by ultra-violet light 
absorption causes atomic oxygen to become increasingly the 
main gaseous constituent. At an altitude of 500 kilometers, 
helium replaces atomic oxygen as the main atmospheric gas. 
At altitudes above 1500 kilometers atomic hydrogen becomes 
the dominant gas (see Figure 7) . 

Traditional research has assigned a drag coefficient of 
2.20 to altitudes between 200 and 400 kilometers. As the 
altitude of the satellite increases above 500 kilometers, 
corresponding to satellite movement through helium, a drag 
coefficient of at least 3.00 must be assigned to maintain 
agreement with observation (see Figure 7) . 

Jacchia [Ref. 7] lists the observed variational changes 
in density and composition in the thermosphere. They are 
the two component solar activity related, semiannual, 
diurnal, seasonal latitudinal, and geomagnetic. 

The two component solar thermospheric variations are 
consequences of changes in solar ultra-violet (EUV) flux. 
The EUV flux changes periodically due to two solar cycles. 
These are the 27 day rotational or disk cycle and the 11 
year sunspot or active area cycle. 

EUV flux is not measurable from the ground due to almost 
complete absorption in the atmosphere. Empirical evidence 
has suggested a strong parallel of the EUV flux and the 
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ground measurable 10.7 centimeter flux (F10.7) which has 
been used in some successful thermospheric models. The 
Jacchia 1960 model used the 10.7 centimeter flux as an input 
to model thermospheric density resulting in an improvement 
of decay predictions over the exponential model. (Indeed, 
it is only recently that there has been a requirement for 
greater accuracy in drag predictions.) 

The thermosphere responds to increased EUV flux by 
absorbing more energy, causing increased heating. This 
increased heating causes an increase in temperature and 
density. The thermospheric response to changes in the EUV 
flux is nonlinear. Increases in EUV flux due to disk 
(rotational) effects produce more thermospheric heating then 
a like increase in flux due to active area effects. 
However, Jacchia noticed in 1973 that the disk component of 
the EUV flux can be linearly related to an average 10.7 
centimeter flux. He recommends using a 71 day average 
corresponding to three solar rotations [Ref. 7]. Tables for 
the 10.7 centimeter flux have been developed.^ The Jacchia 
1977 atmosphere incorporates the average 10.7 centimeter 
flux as an input, used as an indicator of the thermospheric 
heating capacity of the sun's radiation. 



^Orr [Ref. 8: Appendix B] lists values compiled and 
predicted by the NASA Marshall Space Flight Center. Values 
range from a minimum of 73.3 to a maximum of 242.9. The 
units for the 10.7 centimeter flux are 10”22 
watts/meter^/Hz . 
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The semiannual temperature and density variation has 
been modeled in many ways. Jacchia [Ref. 7:pp. 47-50] fits 
a modified harmonic empirical model to the biannual, 
sinusoidal cycle of observed density. Maxima occur in April 
and October. The Jacchia 1971 and 1977 models for the 
semiannual density variation are based on observations of 12 
years of satellite data. The observed density variations 
are thought to be a consequence of the earth-sun 
orientation. The mechanism for the observed effect is 
unknown, but may be related to the maximum equatorial 
heating at the equinoxes. Equinoxes (sun at maximum 
altitude over the earth's equator) occur in late March and 
September. Figure 8 illustrates the semiannual bulge in 
temperature and density observed just after the equinoxes. 
Due to the earth's oblation, the bulk mass of the atmosphere 
is concentrated at the equator. The most direct exposure of 
the equatorial atmosphere occurs at the equinoxes. The 
temperature and density thus peak shortly after periods of 
maximum heating. 

Seasonal-latitudinal thermospheric variations are again 
earth-sun orientation related. Variations in density and 
composition are dependent on the season and latitude. The 
Jacchia 1977 model fits modified harmonic empirical models 
to observed data. The seasonal latitudinal variations 
interact strongly with the semiannual variation, complicat- 
ing the modeling. The winter helium bulge is the best known 
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Figure 8. The Semiannual Atmospheric Density Bulge 

effect of this variational factor. Figure 9 shows the 
general shape change of the temperature and density bulge 
with season. Note that during the June' solstice the bulge 
is centered at about 45 degrees north latitude whereas 
during the December solstice the bulge is centered at about 
45 degrees south latitude. 
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Figure 9. The Seasonal-Latitudinal Atmospheric 
Density Bulge 

The diurnal bulge is shown in Figure 10. It is caused 
by the daytime heating of the side of the atmosphere facing 
the sun. Lagging the actual (most direct) solar exposure by 
about four hours, the heated side shows a marked increase in 
temperature and density. Many current models track and 
predict the position of the diurnal bulge. 
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Figure 10. The Diurnal Atmospheric Density Bulge 



Geomagnetic heating- of the thermosphere is not well 
understood. At times of high geomagnetic activity 
(corresponding to an index of 7) large heating effects are 
produced. The combined effects of geomagnetic and EUV 
heating can produce changes in density of three orders of 
magnitude at 600 kilometers (Jacchia [Ref. 5:p. 82]). 
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In summary, the thermosphere is a very dynamic region of 
the atmosphere. Many of the observed changes require 
further research to complete analytical models. The current 
modified exponential empirical models can be orders of 
magnitude in error in their density predictions. At higher 
earth orbits orders of magnitude error in density is of 
little consequence.^ 

The model chosen for density and compositional effects 
in the atmosphere depends on the importance of the short- 
term cyclic effects. Computations for orbits of months in 
duration and less than 600 kilometers in altitude will be 
subject to large errors using the simple empirical models 
(Jacchia 1966) . For these short missions the more accurate 
Jacchia 1977 model seems more appropriate. (Performance of 
the Jacchia 1977 empirical model has been remarkable with 
density errors rarely exceeding two percent.) Missions 
lasting for decades may be able to use the computationally 
more efficient Jacchia J60 model to predict very long term 
drag effects. 

The one drag force dependent factor not yet discussed is 
the satellite-to-atmosphere relative velocity. This 
relative motion is influenced by three gross factors; local 



^Above 800 kilometers the drag term in the orbital 
decay model used by the Satellite Control Facility in 
Sunnyvale, California is very small. Its main function is a 
catch-all for errors in other perturbating terms. 
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winds, satellite tangential velocity, and atmospheric 
rotation. 

Thermospheric wind has not been measured. Mesospheric 
wind has been measured based on the observation on meteor 
trails. Fleagle [Ref. 4;p. 82] reports meospheric winds 
exceeding 0.150 kilometers per second. This may sound 
large, but must be compared with the satellite's velocity 
(which is on the order of 10 kilometers per second) . 

Satellite tangential velocity is a major factor in the 
satellite-to-atmosphere relative velocity. When combined 
with the rotating atmosphere it accounts for the bulk of 
relative motion. 

Most earth orbit calculations simplify computation of 
perturbation- induced changes in the orbital elements by 
assuming a fixed earth. When considering drag a non- 
rotating earth assumption can lead to a severe error. The 
earth is not fixed. The atmosphere is pulled along with its 
rotation, and lags the earth rotation by a small fraction. 
Figure 11 shows the rotating atmosphere effect. For low 
inclination (prograde) orbits the atmosphere is actually 
following the satellite. In high inclination (retrograde) 
orbits the atmosphere is moving against the satellite. When 
comparing a typical 16 revolution per day satellite 
tangential velocity in a prograde atmosphere to the same 
satellite in a retrograde atmosphere, a 15 percent variation 
in relative velocity can occur. This is compared to the 
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Figure 11. Orbital Inclinations in a 
Rotating Atmosphere 

unknown (but probably less than two percent effect) of 
atmospheric winds. 

A model used by Parks [Ref. l:p. 4] to account for the 
rotating atmosphere separates the relative velocity of the 
satellite in relation to atmosphere into two parts. The 
model is given by: 
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( 11 ) 



Here Vj-q^ is the satellite to atmosphere relative velocity, 
v-j- is the satellite tangential velocity, and 6 is a factor 
accounting for orbital orientation (inclination) with 
respect to a rotating atmosphere (see Figure 11) . 

Equation 11 in Parks [Ref. l:p. 4] models the 

inclination factor in the following way; 

r 2 

6 = [1 - AWg cos i] (12) 

In this model, 6 is the rotation factor, rp the radius at 
perigee, Vp the satellite velocity at perigee, the 

earth's angular rotation rate (assumed constant) , i the 
orbital inclination, and A the ratio of rates of rotation of 
the atmosphere to the earth's rotation. 

The model for 6 used in the propellant longevity model 
neglects the effect of the atmospheric wind, but this should 
introduce only a small error. Recent atmospheric work, 
Jacchia [Ref. 7], suggests a significant gravity wave 
induced motion in the upper atmosphere. As more experimen- 
tal information on this bulk motion becomes available, the 
incorporation of the gravity wave wind may become 
increasingly important. 

The parameter A used by the propellant longevity model 
in the satellite velocity multiplier submodel is probably 
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not constant. Atmospheric theory suggests it may be 
dependent on both altitude and latitude. The range of 
values for the A parameter is assumed to be from a maximum 
of about 1.0 (at low altitudes) to a minimum of around 0.9 
(at high altitudes) . 

The coding of the theory in NSWC TR 83-243 for satellite 
fuel mass decrement did not compute 5 internally. It was 
user supplied input in the original version. As part of the 
mass decrement model validation at the Naval Postgraduate 
School, code internal computation of 6 was added. This 
allows 6 to automatically adjust the satellite-to-atmosphere 
relative velocity Vj-^]^ for the results predicted by the 
theory as perturbation-induced changes to the inclination 
occur. For a long mission life that exhibits significant 
inclination drift, the internal computation of should add 
accuracy to the mass decrement submodel. In place of the 
original user input for 6 , the A is now specified by an 
input file. The input file default value for A is 0.9. 

A model atmosphere used for drag computations must be 
able to correctly predict the effects on density, drag 
coefficient, and satellite-to-atmospheric velocity. There 
are many variations in the atmosphere that affect density. 
The sum of all these effects creates a density bulge in the 
atmosphere, the location and shape of which follow maximum 
heating. 
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The drag coefficient is atmosphere composition 
dependent. Thus the model atmosphere must be able to 
predict the dominant species in the satellite environment. 

Satellite-to-atmosphere relative velocity is a function 
of satellite tangential speed, atmospheric rotation, and 
orbital inclination. 

The Jacchia J60 model atmosphere used in the propellant 
longevity model only includes the effects of diurnal heating 
on the density bulge. There are no computations of 
atmospheric composition. The newer Jacchia J77 model 
incorporates all the known density variations currently 
observed. The Jacchia J77 atmosphere also predicts 
thermospheric composition. 

The background concepts for the propellant longevity 
model understanding and operation have been presented. The 
next chapter describes the model and its computer program. 
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III. PROPELLANT LONGEVITY MODEL DESCRIPTION 



This chapter describes the propellant longevity model. 
First, a brief summary of the model theory is presented. The 
second part of this chapter describes the coding of the 
model as installed at the Naval Postgraduate School. 
Finally, the third section gives the model code compiling 
and handling instructions specific to the system at the 
Naval Postgraduate School. 

A. PROPELLANT LONGEVITY MODEL THEORY 

A detailed development of the propellant longevity model 
is presented by Dr. Parks [Ref. 1] . The propellant 
longevity model hinges on the prediction of the fuel mass 
used by a satellite to compensate for atmospheric drag. The 
mass decrement submodel predicts the amount of fuel consumed 
after each revolution of the satellite. 

The mass decrement submodel for a low-earth-orbit 
satellite in an oblate diurnal atmosphere is presented in 
Reference 1 as Equation 56. The diurnal density bulge is 
the only density and compositional atmospheric variation 
considered in the model. 

The integrals of the mass decrement model (Equation 56) 
are next expressed in terms of Bessel functions of the first 
kind and order n as defined by: 
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In(eae) = ^ / cos(nE) exp(3ae oosE) dE ( 13 ) 

JIT 0 

Here is the Bessel fxonction of the first kind and order 
n, 3 the inverse density scale height (assumed constant) , a 
the semi-major orbital axis, e the orbital eccentricity, n = 
0,1,2, ..., 6, and E the orbital eccentric anomaly. All of 

these orbital mechanics terms were presented in Chapter 
I I. A. The integral is to be evaluated over one revolution. 

The Bessel function expansion of the mass decrement 
model is given in equation 58 of Reference 1. This 
particular expansion is central to the model. Provided 
accurate values for other multipliers are produced by the 
model atmosphere, (accurate values for atmospheric density, 
satellite drag coefficient, and satellite-to-atmosphere 
relative velocity are essential) greatly increased accura- 
cies should result in predicting the fuel consumed to 
compensate for atmospheric drag. The main draw back to the 
use of the this type of Bessel function is the known 
slowness of its convergence. 

B. COMPUTER CODE DESCRIPTION 

A detailed description of the model code and its input 
elements is contained in the code and included in Appendix 
A. An operational description is presented here. 

The computer code calculates the propellant mass 
remaining together with the cumulative mission life after 
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each revolution of the satellite. When propellant fuel is 
exhausted, the computation terminates. 

Computer output is available in several forms and is 
specified by the user through the use of some control 
parameters. The user can select one of two modes; standard 
mode or a sensitivity analysis mode. The input elements used 
as control parameters for output and mode control are: 
maximum number of revolutions ITERl, sensitivity analysis 
granularity PSIZE, output print frequency ITER2, mode 
control PSENS, anomalistic period print control lET, number 
of scheduled orbit adjusts NOA, scheduled orbit adjust array 
IRV(K) , and the eimovint of each adjust array DA(K) . These 
terms are discussed below. 

The most basic form of output available is a tabulation 
of propellant remaining and accumulated mission life for 20 
specified satellite-orbital initial conditions. The 20 
initial conditions are specified in an input file. 
Accumulated mission life is expressed in modified Julian 
days (MJD)^ and revolution number. Frequency of output 
printing is selectable from a maximum of once per revolution 
to a minimum of once per 10000 revolutions by user 
specification of ITER2 . 

A cap on runaway computer time is available by selection 
of ITERl. Satellites with small cross-sectional areas, high 

^A description of the modified Julian day time 
measurement is in Taff [Ref. 2:p. 103]. The modified Julian 
day is reckoned from 0^ Nov. 17 1858. 
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altitudes, and large amounts of maneuvering fuel can have 
very long orbital life-times. By defining a maximum number 
of revolutions to be considered in the calculations through 
the reasonable selection of ITERl, unnecessary expenditure 
of computer time can be avoided. 

As part of the model verification a sensitivity analysis 
option was added to the program. It can be selected by 
setting PSENS = 1. The standard operating mode is selected 
if PSENS =0. In the sensitivity analysis mode, the 
subroutine SNALYS creates an array of elements selected by 
the user to be linearly scanned from a specified minimum to 
maximum with a given granularity. 

Granularity is specified by the selection of the control 
parameter PSIZE. A granularity of 1/100 (code specified as 
PSIZE = 100) will cause 100 executions of the propellant 
longevity code to complete fuel exhaustion. This 
specification can lead to excessive computer execution times 
for high altitude satellites. A granularity of 1/5 will 
cause five executions of the code to fuel exhaustion. 
However, only five iterations from a given minimum to a 
given maximum may be too coarse to reveal the true trend of 
the satellite mission life. 

Each time the code is sequentially executed in the 
sensitivity analysis mode, the elements selected by the user 
for analysis are incremented linearly. The output in the 
sensitivity analysis mode is a table of mission lifetimes. 
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resulting from the initial values of input elements not 
selected for analysis, and the current value of the elements 
selected for analysis. 

The output always echoes the initial input element 
conditions as read from a user supplied input file. Further 
choice of output and mode control parameters allow the user 
to tailor the output. 

Three more input control elements allow for the 
specification by revolution number of up to ten adjusts to 
the semi-major axis. The element IRV(K) is an array whose 
elements specify the orbital revolution number for axis 
adjust. The element DA(K) specifies the amount of change in 
the semi-major axis in kilometers for each scheduled 
adjustment. Finally, NOA specifies the total number of 
orbital adjusts to be performed. If the user does not 
desire any changes to the semi-major axis, NOA must be set 
to 0 and no values for IRV(K) and DA(K) input to the file 
stream. 

An additional twenty input elements must be specified by 
the user in the input file. Five of these are satellite 
physical parameters, five are atmosphere specific 
parameters, and ten are the Brouwer mean orbital element set 
of the NAVSPASUR system. 

The five satellite physical parameters are: propulsion 
motor specific impulse (SPI in seconds) , initial mass of the 
satellite including fuel (WT in kilograms) , initial mass of 
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the maneuvering fuel (W in kilograms) , satellite atmospheric 
drag coefficient (CD) , and satellite cross-sectional area (A 
in square kilometers) . 

Motor specific impulse SPI can range from 220 seconds 
for monopropellants to 3000 seconds for ion drives. The 
default value is 230. 

The initial satellite mass WT depends on the satellite 
mission and design. A typical communications satellite has 
a mass of about 3500 kilograms. An SDI satellite may be 
considerably more massive. The default value for the 
satellite mass is 9100 kilograms. 

Initial maneuvering fuel mass W of the satellite depends 
on the design life of the satellite as well as it's 
ballistic coefficient. Low altitude satellites need at 
least one percent of their initial mass as maneuvering fuel 
if mission lives of years are to be achieved. The default 
value for initial maneuvering fuel is 100 kilograms. 

The drag coefficient CD is a measure of the slowing 
power of the atmosphere on a given satellite. In the 
thermosphere it is very atmosphere composition dependent. A 
value of 2.2 has been settled on by orbital analysts for 
altitudes of 200 to 400 kilometers. A value of at least 3.0 
must be assigned at altitudes above 500 kilometers 
(corresponding to the helium atmosphere) . The default value 
for the drag coefficient is 2.2. 
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Satellite cross-sectional area is again very design 
dependent. The echo series satellites had cross-sections 
measuring in hundred's of square meters. Some SDI 
satellites may also have very large cross-sections. The 
default for cross-sectional area is 50.0*10”^ square 
kilometers, corresponding to 50 square meters. 

The five atmosphere specific parameters are: atmosphere- 
to-earth relative rotation rate D, 10.7 cm solar flux F107 
described in Chapter II, ninety-day average solar flux FBAR, 
geomagnetic activity index AKP, and the time of vernal 
equinox passage TVE. 

The atmosphere-to-earth relative rotation rate D is the 
factor explained at the end of Chapter II. The default 
value is 0.9. 

The decimetric solar flux and the average decimetric 
solar flux F107 and FBAR are again explained in Chapter II. 
Default values of 100 are assigned to each. 

The geomagnetic index AKP is discussed in Chapter II. 
It ranges from 0 to 7 and is set at a default value of 2.00. 

With the current atmospheric model Jacchia J60, both 
FBAR and AKP are read as input but not used for the density 
computation. As seen from a sensitivity analysis, these 
parameters currently have no effect on mission life. The 
Jacchia J77 model atmosphere uses these now stranded input 
parameters. The future incorporation of the Jacchia J77 
atmosphere model will be able to use these input elements. 
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The ten NAVSPASUR orbital elements are the final input 
elements in the standard input file. They are all explained 
in the orbital mechanics section, and in the code included 
in Appendix A. Briefly, they are; 



UJD; 


time of epoch in modified julian 


days 


ETU; 


initial mean anomaly 




GO; 


argument of perigee 




HO; 


right ascension of the ascending 


node 


ES; 


eccentricity 




XI; 


inclination 




AO; 


semi-major axis 




RND; 


rate of change of mean motion 




ESD; 


eccentricity decay rate 




ADOT; 


semi -major axis decay rate. 





The first seven of these are orbital elements and the final 
three are the results of perturbations. The NAVSPASUR 
orbital element set is the result of orbital measurements of 
actual satel 1 ites . 



The final block 


of input in the 


input 


file 


is 


the 


sensitivity analysis 


specification array. 


It 


must 


be 


included only when 


the sensitivity 


analysis 


mode 


is 



specified. The array consists of 20 lines used to make up 
the array of elements to be scanned by the program when the 
sensitivity analysis mode is selected. 

The sensitivity analysis specification array is divided 
into four columns. The first column contains the input 
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element variable name. These are the same names as 
specified in the satellite physical, atmosphere specific, 
and orbital elements input discussions. 

The second column is a control parameter. When set to 
1.0, the corresponding input element is included in the 
sensitivity analysis. Colximns three and four are in a 
Fortran F16.9 format and they specify the minimum and 
maximum range of the element to be scanned. Defaults for 
the minimum and maximum have been provided in specific cases 
corresponding to logical and physical ranges of the input 
elements. In some cases zeros appear in the minimum' and 
maximxim columns. These zeros correspond to input elements 
not physically constrained. An example of an input element 
not physically constrained is the time of epoch (UJD) . 

If a 1.0 is entered in column two of a particular row of 
the sensitivity analysis specification array, that input 
element is included in the analysis. Any number of the 
input elements may be selected for simultaneous analysis. 
When selecting a multivariable sensitivity analysis, the 
user must understand the interactions of the incrementing 
variables and ensure that a solution exists in the mission 
life range being considered. 

When a variable is included in the sensitivity analysis, 
the incrementing values override the initial conditions. 
All input elements not user selected for analysis (those 
with zeros in the second column) are reset to their initial 



60 



values after each iteration during the analysis. Appendix B 
contains sample input element files and their corresponding 
outputs . 

Input elements and output forms have been explained. 
The program flow is discussed next. The 16 subroutines and 
eight functions comprising the propellant longevity model 
are explained in Appendix A in the Fortran code. The 
routine PLEP initializes the system, including inputting 
program control parameters and initial values for all 
elements. If the sensitivity analysis mode is selected, 
subroutine SNALYS generates the array to be scanned 
linearly, complete with the necessary increments for each 
iteration. Analysis mode causes completion of all 
computations to propellant exhaustion PSIZE times. The 
routines SNALYS and PLEP are completed only once even in 
analysis mode. 

Program control is transferred from the initialization 
routine PLEP to the subroutine KOZAI, the main driver and 
bookkeeper of the system. The routine KOZAI first causes 
the NAVSPASUR Brouwer mean element set to be Converted into 
a Keplerian (Kozai-like) mean orbital set, by calling the 
subroutines BRAUER, PERIOD, and MEAN. 

S\abroutine MEAN calls functions XSPA-XSPM to convert the 
Brouwer osculating elements generated in BRAUER to Keplerian 
orbital elements. The Keplerian elements are then returned 
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to the subroutine KOZAI. In standard mode the mean element 



set is given as output. 

All of the above computations have simply conditioned 
the input elements to conform to the requirements of the 
propellant longevity model. The routine KOZAI now causes 
calculation of the fuel used to compensate for atmospheric 
drag after each revolution. The computation is continued 
until the fuel is exhausted or the maximum number of 
revolutions (ITERl) is reached. When the computation 
terminates, program control is transferred back to PLEP. 

To calculate the drag-compensating fuel required for 
each revolution, KOZAI calls the subroutines PERIOD, GEOP, 
THRST, and DRAG. 

The subroutine GEOP calculates the change in orbital 
elements due to the aspherical earth. It uses the J2, J3, 
and J4 terms in a Legendre polynomial expansion, as 
explained in Chapter II. Aspherical earth calculations are 
performed once each revolution. 

The subroutine THRST checks each revolution for a 
scheduled orbital adjust. If a change to the semi-major 
axis is scheduled, it calculates the change to the axis, the 
argument of perigee, the eccentricity, and the mean anomaly. 
It also computes the amount of maneuvering fuel required to 
perform the adjustment. All changes are returned to the 
subroutine KOZAI. 
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The siibroutine DRAG calculates the fuel mass needed to 
overcome the atmospheric drag. This calculation is 
accomplished by calling four subroutines: SOLOC, SATLOC, 
SALT, and JAC60. 

The routine SOLOC computes the right ascension and 
declination of the sun for the model atmosphere diurnal 
density bulge. The routine SATLOC computes the right 
ascension, declination and geomagnetic latitude of the 
satellite. This is accomplished by first calling the 
svibroutine POSVEL to calculate satellite position and 
velocity. The routine POSVEL in turn calls the subroutine 
NWTRPH to solve Kepler's equation (see Chapter II). The 
routine SALT computes the satellite altitude. 

Atmospheric density is computed by the subroutine JAC60 
based on the Jacchia J60 model atmosphere discussed in 
Chapter II. A maximiom, minimum, and mean density that the 
satellite will encounter in the oblate diurnal atmosphere 
are calculated. 

The mass decrement equation (Equation 58 of Reference 1) 
is contained in the subroutine DRAG. The function MMBSIR is 
used to perform the Bessel function expansion. The routine 
DRAG returns the fuel mass decrement due to atmospheric drag 
computation to KOZAI. The routine KOZAI then decrements the 
remaining fuel by this amount and checks for fuel exhaustion 
and maximum orbital number. If fuel is less then zero, or 
the current revolution number equals ITERl, calculations are 
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terminated and control is returned to PLEP. 



In standard 



mode, program execution now halts. 

In sensitivity analysis mode, all input elements not 
selected for analysis are returned to their initial values, 
elements selected for analysis are linearly incremented, and 
KOZAI is again called by PLEP. The subroutine KOZAI then is 
executed, as previously described, with the new values for 
the selected elements. Upon fuel exhaustion or reaching the 
maximum revolution number, KOZAI terminates calculation and 
passes program control back to PLEP. This process continues 
until the elements user selected for analysis have been 
incremented from their minimum to maximum. When the maximum 
value of the elements user selected for analysis is reached, 
program execution terminates in the sensitivity analysis 
mode. 

The tabular output in the analysis mode is printed line- 
by-line in KOZAI upon reaching fuel exhaustion or maximum 
revolution number. The result is a table of data, PSIZE 
long, that represents the linear scan of the specified 
elements. The first line is the mission life corresponding 
to the minimum value of the scanned elements. Each 
subsequent line represents the total mission life resulting 
from an increment of the elements selected for analysis and 
the initial conditions of the elements not selected. The 
final line is the mission life resulting from selected 
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elements at their maximum value and non-selected elements at 
their initial conditions. 

Mission life is expressed in revolution number and 
modified Julian days. The last column of the tabular output 
in analysis mode represents the current value of the 
incremented input element. When more than one element is 
analyzed, the last column is the current fraction of scan, a 
real number between 0 and 1. 

The tabular output can easily be converted into an input 
file for a graphics routine by deleting the first 98 lines 
of output. 

C. MODEL CODE OPERATION 

This section describes the operation of the computer 
code at the Naval Postgraduate School on the IBM 3033 
computer under VMS using the VS FORTRAN compiler. 

The propellant longevity model must be compiled in the 
VM mode using the automatic precision increase option. If 
the propellant longevity code is in filename SDI and 
filetype FORTRAN (the filetype must be FORTRAN) , the correct 
form for the compilation command is; 

FORTVS SDI (AUT0DBL(DBLPAD4) ) 

The listing generated by this compilation takes up about 20 
percent of the user's "A" disk. 

Prior to program execution, the input and output files 
must be specified. If the input file is in filename PLEPOl 
and filetype DATA2, then the correct form of the command is: 
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FILEDEF 5 DISK PLEPOl DATA2 A1 



This informs the computer where to look for input during 
execution of the model code. If hard copy of the output is 
desired, it can be written to a file instead of the screen 
by using the FILEDEF command. If the file name is SDI and 
the desired type is DATA, the correct command is: 

FILEDEF 6 DISK SDI DATA A1 

After execution of the propellant longevity code, the 
results will be in SDI DATA on the user's "A" disk. 

To execute the model listing generated upon compilation, 
the command RUN SDI is typed (this assumes that SDI is the 
filename of the propellant longevity program) . 

This program can have very long execution times. A 
general rule of thumb is that execution time takes less than 
PSIZE * ITERl/10,000 minutes. Roughly, the computation 
takes 100 microseconds for each revolution. For example, if 
PSIZE is 100 and ITERl is 10,000 (corresponding to a five 
year lifetime) then the IBM 3033 will take 200 minutes to 
complete the calculations. The user must carefully choose 
small values of PSIZE and ITERl on the initial investigation 
of elements. 

The simplest graphical routine to construct plots of 
tabular data is EASYPLOT. The EASYPLOT routine was used to 
generate all the plots shown in Chapter IV. 

In summary f the steps of program operation are: 

(1) Compile the program using the automatic 
precision increase option. 
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(2) Select the mode of operation and input element 
values by manipulation of the input file. 

(3) Define input and output files. 

(4) Execute the program by typing the file name. 
The next chapter presents the results of the 

element sensitivity analysis performed to test the 
reasonableness during verification. 



input 

model 
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IV. MODEL VERIFICATION 



This chapter presents the results of the propellant 
longevity model verification at the Naval Postgraduate 
School. Standard model verification tests the model against 
three questions. First, does the model address the problem? 
Second, are the model predictions reasonable? Third, how do 
the model predictions compare to real-world observations? 

This verification addresses only the first two criteria. 
Observational data for comparison with the model predictions 
is not available in an unclassified state. The quantitative 
analysis of model predictions is thus not addressed in this 
verification. 

The model does address the initial problem. It predicts 
the fuel mass decrement of a low-earth- orb it satellite doing 
intrack micro-thrusting to overcome atmospheric drag. 

The reasonableness test is the subject of the remainder 
of this chapter. An effective way to test the model 
performance is to examine model predictions resulting from 
the full ” intended use” range of the input variables. 
Furthermore, if a plot is made of predicted results as a 
selected input element is incremented through its range, 
prediction trends can be found and examined. This procedure 
provides a qualitative analysis of model sensitivity to the 
selected input .element. 
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To aid in these tests and trend studies, a sensitivity 
analysis feature was incorporated into the program code of 
the model, as previously described. Graphical results of 
each single variable sensitivity analysis are presented for 
the input elements . 

The format of the following figures are all related. 
They represent plots of satellite mission life in modified 
Julian days (here after referred to simply as "days”) versus 
a specified input element. All figures have the Y-axis as 
mission life and the X-axis as the variable selected for 
analysis. In each case the range of values for the input 
element is selected to remain within reasonable values of 
the input variable, as discussed in Chapter III. The 
granularity in most cases is 1/25 (PSIZE = 25) . In some 
cases it is necessary to decrease the granularity to 1/100 
(PSIZE = 100) to see the actual trends. 

The default values discussed in Chapter III are used for 
all input elements not selected for analysis (except where 
noted) . They correspond to a near-polar, retrograde orbit 
of an actual satellite supplied as sample input by NSWC. 
The model predicted propellant longevity using all default 
conditions and 10 kilograms of fuel is 13.74 days. When the 
default values and 100 kilograms of fuel are used, the 
resulting propellant longevity is 150.21 days. An * symbol 
on the plotted curves indicates the default value of the 
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analyzed input element. The default value names are again 
summarized here. The units can be found in Chapter III. 



Satellite physical parameters; 
SPI (motor specific impulse) 

WT (satellite initial mass) 

W (initial fuel mass) 

CD (drag coefficient) 

A (cross-sectional area) 
Atmospheric specific parameters; 
D (atm. rotation factor) 

F107 (decimetric solar flux) 
FBAR (average F107) 

AKP (geomagnetic index) 

TVE (time of vernal equinox) 
NAVSPASUR orbital elements; 

UJD (time of epoch) 

ETU (mean anomaly) 

HO (R.A. of ascending node) 

GO (argument of perigee) 

ES (eccentricity) 

XI (inclination) 

RND (time rate of mean motion) 
ESD (eccentricity decay rate) 
AO (semi-major axis) 

ADOT (semi-major axis decay) 
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The input element sensitivity analyses are divided into 
related sets. All the sets are discussed in the same basic 
format. Each analysis is single-variable and arranged as 
follows: 

(1) The expected trend of real-world mission life 
resulting from the selected variable 
incrementing through its range. 

(2) The expected mission life trend predicted by 
the model as a result of the selected 
variable incrementing through its range. 

(3) The actual mission life results of the 
selected variable incrementing through its 
range . 

(4) A discussion of the selected variable 
sensitivity analysis. 

The explanations for the plotted trends are based on 
expected behavior in agreement with current theory. An 
understanding of the actual mechanisms producing the 
observed trends will be improved with further research. 

For a given set of satellite-orbital configurations, the 
satellite mission life will be shortest for satellites 
traveling through the most dense atmosphere. Thus for a 
given input variable sensitivity analysis, the range of 
values corresponding to satellite travel through the most 
dense atmospheric conditions will produce the shortest 
mission life. The range of values of the same input 
variable corresponding to satellite travel through the least 
dense atmosphere will likewise produce the longest mission 
life. The orbital-perigee atmospheric-density-bulge 
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encounter geometry is thus critical in mission life 
analysis. 

The initial sensitivity plot of an input variable is 
accomplished with 10 kilograms of fuel. In many cases, true 
geopotential-induced orbital-drift trends are not seen with 
small amounts of drag-compensating fuel. To verify trends 
caused by geopotential-induced orbital drift, the initial 
compensating fuel is increased to 100 kilograms. The 
differences in the 10 kilogram and 100 kilogram mission life 
plot trends can then be attributed to geopotential effects 
and the movement of the center of most direct solar heating 
during the mission. 

The 133.14 day nominal mission life resulting from the 
default values of the input elements and 100 kilograms of 
drag compensating fuel is almost a half-year. During this 
half-year, the earth is revolving around the sun causing the 
noon-sun atmospheric position to increase in longitude by 
about one degree per day. This greatly complicates trend 
analysis when using the 100 kilogram fuel case alone. It is 
thus necessary to include the 10 kilogram fuel case. The 
13.74 nominal mission life of the 10 kilogram fuel case is 
short enough so that it shows little longitudinal movement 
of the most-direct-sun atmospheric position. 

In many of the sensitivity analysis plots, the predicted 
mission life curve has a "knee." This "knee" corresponds 
to the range of the selected input variable that produces 
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changing behavior in the mission life curve. The location 
of these ranges has significance to design research. 

The input variables are divided into four categories 
based upon their effect on mission life. The first category 
contains the input variables that affect the atmospheric- 
bulge orbital-perigee encounter geometry. The second 
category contains the input variables that affect the 
slowing effectiveness of the atmospheric density bulge on 
the particular satellite configuration. The third category 
contains the input variables that affect the depth of 
penetration of the satellite into the atmosphere at perigee. 
Finally, the fourth category contains the input variables 
that have no affect at all on mission life. As will be 
seen, predicted mission life is most sensitive to the third 
category (factors affecting atmospheric penetration at 
perigee) . 

Factors affecting the atmospheric-density-bulge orbital- 
perigee encounter geometry are: the time of vernal equinox 
passage TVE, the time of epoch UJD, the argument of perigee 
GO, the right ascension of the ascending node HO, and the 
inclination XI. The time of vernal equinox provides a real 
time reference for earth in its orbit relative to the sun. 
When combined with the time of epoch (sometimes called the 
time of perigee passage) , these two factors determine the 
position of the noon sun in right ascension and declination 
relative to the atmosphere. The default values used in the 
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model verification position the initial noon sun at 
approximately 300 degrees longitude and 20 degrees south 
latitude relative to the celestial meridian containing the 
first point of aries. The right ascension of the ascending 
node and the argument of perigee combine to position the 
point of perigee in longitude relative to the first point of 
aries. The orbital inclination positions the point of 
perigee in latitude (or declination) . Actual positioning of 
the satellite in the above orbit is accomplished by 
specifying the mean anomaly. All these orbital mechanics 
terms have been explained in Chapter I I. A. 

The expected real-world behavior of satellite mission 
life as a result of varying the time of vernal equinox or 
time of epoch is discussed next. The time of vernal equinox 
TVE serves to set a real-time reference for positioning of 
the earth in its orbit relative to the sun. When the time 
of epoch UJD is specified, the position of the point of the 
most direct atmospheric heating (noon sun) in the atmosphere 
is set. The atmospheric density-bulge follows the location 
of the most direct solar heating. The result of TVE and UJD 
varying in relation to each other is a seasonal positioning 
of the earth in its orbit around the sun at the time of 
initial perigee passage. Due to the earth's rotational-axis 
tilt, the point of most direct solar heating moves in 
longitude and latitude as the earth moves in its orbit 
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around the sun. Figures 3 and 4 in Chapter II. A aid this 
visualization. 

The earth and solar cycles effecting the atmospheric 
density-bulge combine to define its size and shape. These 
cycles are discussed in Chapter II. B, and stated again here. 
They are: diurnal, semi-annual, seasonal-latitudinal, solar 
rotational, and solar activity related. The Jacchia J60 
model atmosphere, as used in the propellant longevity model, 
incorporates only the diurnal cycle. 

Once the position, size, and shape of the atmosphere is 
set, the mission life of a drag compensating satellite 
orbiting within it can be developed. For such a satellite 
in a specific initial orbit, the following real-world trend 
in mission life (propellant longevity) should be observed. 
The trend should be cyclic, resulting from the superposition 
of all the atmospheric heating cycles. Mission lives will 
be shortest for orbits having their perigee in the most 
dense atmospheric bulges. The longer mission lives will be 
observed for satellites having their orbits farthest away 
from the density-bulge. Thus the plots of satellite mission 
life will show peaks and valleys consistent with the above 
discussion. 

In the propellant longevity model, the diurnal cycle is 
the only atmospheric heating cycle that is incorporated. 
Positioning of the earth in its orbit is accomplished by the 
combination of the time of vernal equinox and time of epoch 
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in the subroutine SOLOC. As discussed in Chapter III, the 
subroutine SOLOC positions the sun by right ascension and 
declination over the earth's atmosphere. The routine DRAG 
considers the diurnal density bulge to be centered under the 
solar position. The density profile of the atmospheric 
density-bulge is generated in the routine JAC60 based on the 
value of the input element F107 . 

The propellant longevity model considers the atmosphere 
to be oblate, diurnal, and revolving around the sun with a 
23.5 degree rotational-axis tilt. The expected model- 
predicted-mission-life trend as is cyclic as TVE and UJD are 
varied. Every year should be an exact replica of every 
other year. Two peaks and two valleys in the yearly 
mission-life plot should be observed. The deepest valley 
occurs when the point of perigee is most closely centered in 
the atmospheric density bulge. With a low eccentricity (e = 
.01) orbit considered, a much shallower valley will occur 
about six months later. This shallower valley corresponds 
to the point of apogee being most . closely centered in the 
atmospheric density-bulge. 

Two peaks in plotted mission life should occur 
corresponding to the orbital pane being most perpendicular 
to the earth-sun line. The peaks should be six months 
dpart, of approximately equal magnitude, and interspaced 
between the mission-life valleys. As the eccentricity of 
the orbit is increased, the apogee-associated valley will 
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become shallower while the perigee-associated valley becomes 
deeper. 

Figures 12, 13, 14, and 15 present the plots of the 
sensitivity analysis of TVE. In Figure 12, TVE has a range 
of 1000 days, from 38422 to 39422 days. The initial fuel is 
10 kg. The resultant predicted mission life ranges 
cyclically between 11 and 18 days. 

Figure 13 has an incremented range for TVE of 38422 to 
38522, or 100 days. Initial fuel is 10 kg. The resulting 
predicted mission life again ranges cyclically between 11 
and 18 days. 

In Figure 14 the width of the incremented range is 
raised to 10,000 days (about 27 years). The minimum value 
of TVE is 38422 and the maximum is 49422. Again 10 kg of 
fuel is used. The results are cyclic between 11 and 18 days 
until TVE becomes greater than time of epoch UJD. The 
results at the high end of TVE in this plot are thus 
unreasonable. The time of vernal equinox passage chosen 
must occur before the satellite is launched requiring TVE to 
be less then UJD. 

The TVE is again incremented through a small range of 20 
days in Figure 15. Initial fuel is 100 kg. The resulting 
range of mission life predicted is from 150 to 153 days. 
Note the absence of the cyclic trend*. 

Figures 16 and 17 present the results of the sensitivity 
analysis plots of the time of epoch UJD with a range of 1000 
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Figure 12. Mission Life versus Time of Vernal 
Equinox (1000 Days) 
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MISSION LIFE (DAYS) 




Figure 13. Mission Life versus Time of 
Vernal Equinox (100 Days) 



79 



MISSION LIFE (DAYS) 




Figure 14. Mission Life versus Time of 
Vernal Equinox (10,000 Days) 
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Figure 15. Mission Life versus Time. of 
Vernal Equinox (20 Days) 
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MISSION LIFE (DAYS) 




Figure 16. Mission Life versus Time of Epoch 
(1000 Days and 10 kg of Fuel) 
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Figure 17. Mission Life versus Time of Epoch 
(100 Days and 100 kg of Fuel) 
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days. Note the similarity to the TVE plots, as 
theoretically expected. Figure 16 shows the 10 kg fuel 
case, and Figure 17 the 100 kg fuel case. 

The trends observed in Figures 12 through 17 are in 
agreement with expected model predictions. The plotted 
trends may have significant differences from real-world 
trends due to the limitations of the Jacchia J60 model 
atmosphere. It should be noted that the default values of 
TVE and UJD produce a splar location approximately 30 
degrees longitude past the winter solstice in the southern 
hemisphere heading toward the spring (vernal) equinox. The 
approximate latitude of the default solar position is 20 
degrees south . 

The right ascension of the ascending node HO positions 
the orbital plane with reference to the celestial meridian 
containing the first point of aries. When the argument of 
perigee GO and the inclination XI are specified, the point 
of perigee is positioned in longitude and latitude. This 
geometry is shown on Figure 4 in Chapter II. A. 

The real-world trend in mission life of a drag 
compensating satellite, as a result of varying HO, GO, and 
XI is cyclic. The cycle repeats every 360 degrees for HO 
and GO, and every 180 degrees for XI. The observed pattern 
of peaks and valleys will be highly variable, depending on 
the values of these three factors. 
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Model predicted mission life trends due to varying HO, 
GO, and XI through their ranges should closely parallel 
real-world expectations. This close parallel between model 
predictions and observation assumes the atmospheric density- 
bulge has been properly positioned and described. 

The default values of 95 and 200 degrees for the 
inclination and the argument of perigee combine to produce 
a point of perigee that is very nearly 20 degrees below the 
descending node in the southern hemisphere of the earth's 
atmosphere. The results of the sensitivity analysis plots 
for the right ascension of the ascending node and the 
argument of perigee are shown in Figures 18 through 21. 

The right ascension of the ascending node HO is 
incremented from 0 to 360 degrees. Figure 18 presents the 
10 kg fuel case and Figure 19 the 100 kg fuel case. 
Predicted mission life ranges are 11 to 18 days in Figure 
18, and 135 to 160 days in Figure 19. 

The argument of perigee GO is incremented from 0 to 360 
degrees with a 10 kg fuel case in Figure 20, and a 100 kg 
fuel case in Figure 21. The predicted mission life are 13 
to 18 days, and 148 to 160 days, respectively. 

In Figure 18 two approximately equal mission life peaks 
and valleys are observed. The valleys correspond to the 
points of perigee and apogee closely approaching the 
atmospheric density bulge. The peaks occur when the orbital 
plane is most nearly perpendicular to the earth-sun line 
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Figure 18. Mission Life versus Right Ascension 

of the Ascending Node (10 kg of Fuel) 
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Figure 19. Mission Life versus Right Ascension 

of the Ascending Node (100 kg of Fuel) 



87 





ARGUMENT OF PERIGEE (DEGREES) 



Figure 20. Mission Life versus Argument of 
Perigee (10 kg of Fuel) 
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Figure 21. Mission Life versus Argument 
of Perigee (100 kg of Fuel) 
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which contains the atmospheric density bulge center. The 
change in shape of the 10 and 100 kilogram curves (Figures 
18 and 19) is due to the movement of the atmospheric density 
bulge during mission life and to geopotential-induced 
orbital drift. 

With default values of 95 degrees for both the right 
ascension of the ascending node and inclination, the effect 
of varying the argument of perigee from 0 to 360 degrees is 
to move the point of perigee around the initially fixed 
orbital plane (see Figure 4 in Chapter II) . The plot of 
mission life corresponding to this motion should by cyclic 
showing one peak and one valley. The valley corresponds to 
the perigee point closest to the atmospheric density-bulge 
center. The peak in plotted mission life occurs when the 
apogee point is closest to the atmospheric density-bulge 
center. 

Figures 20 and 21 show the plotted results of the 
sensitivity analysis of the argument of perigee GO. The 
expected single peak and single valley are seen. The 
difference in the 10 kilogram case in Figure 20 and the 100 
kilogram case in Figure 21 is due to the seasonal movement 
of the atmospheric density-bulge and geopotential- induced 
orbital drift. The magnitude of the mission life variations 
is less than that produced by varying the right ascension of 
the ascending node. • This effect is due to low eccentricity 
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orbit (e =0.01) and the position of the orbital plane with 
reference to the earth-sun line. 

Factors which primarily move the orbital plane in 
relation to the atmospheric density-bulge center will show 
two peaks and two valleys in their single-cycle mission-life 
range. Factors which primarily move the point of perigee in 
an initially specified orbital plane will show one peak and 
one valley in their single-cycle range. These 
considerations assume that only low eccentricity orbits are 
being considered. 

Based on the above discussion, Figures 18 through 21 
appear to agree with expected real-world and model- 
predicted mission life trends. 

The analysis of predicted mission life trends due to 
varying the inclination XI is more complex. Several factors 
must be considered in addition to the orbital-perigee 
atmospheric density-bulge encounter geometry. First, the 
diurnal atmosphere model used by the propellant longevity 
model is also oblate, and is thus more dense at the equator 
than at the poles. This fact means that satellites in near 
polar orbits travel through an overall less dense atmosphere 
than those in. a near equatorial orbit. Second, the 
atmosphere rotates in the same direction as the earth. 
Hence prograde orbital satellites experience less satellite- 
to-atmosphere . relative velocity than satellites • in 
retrograde orbits. If the satellite altitude is low enough. 
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a near-equatorial prograde orbit satellite has a longer 
mission life than the same satellite in a near-equatorial 
retrograde orbit. These additional effects complicate the 
analysis of mission life resulting from varying the orbital 
inclination. 

The default value position of the center of the 
atmospheric density-bulge is 300 degrees longitude and 20 
degrees south latitude. The bulge is moving at a rate of 
about one degree longitude per day toward the vernal 
equinox. The point of perigee is default positioned at 
approximately 275 degrees longitude and 20 degrees south 
latitude. As can be seen, the point of perigee is quite 
close to the center of the atmospheric density-bulge. The 
motion of the point of perigee, as the inclination is 
increased from near 0 degrees to near 180 degrees, describes 
a semicircle 20 degrees in radius around the descending 
node. The point of perigee remains quite close to the 
initial position of the center of the atmospheric density- 
bulge . 

Figures 22, 23, 24, and 25 are the predictions of the 
sensitivity analysis plots of orbital inclination XI. The 
incremented range is from 22.5 to 157.5 degrees. Nearer 
equatorial inclinations are prohibited by the available 
arithmetic precision on the IBM 3033. Figure 22 shows the 
10 kg fuel case. Figure 23 the 100 kg fuel case, and Figure 
24 the 100 kg fuel case with a granularity of 1/100 and a 
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Figure 22. Mission Life versus Orbital 
Inclination (10 kg of Fuel) 
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Figure 23. Mission Life versus Orbital 
Inclination (100 kg of Fuel) 
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Figure 24. Mission Life versus Orbital Inclination 
(100 kg of Fuel and PSIZE = 100) 
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Figure 25. Mission Life versus Orbital Inclination 

(100 kg of Fuel, PSIZE = 100, and D = 0.9) 
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rotation factor of 0.9. Predicted mission life ranges are 
7.8 to 13.7 days in Figure 22, and 90 to 138 days in Figures 
23 and 25. 

Figure 25 shows the effect of the using 100 kilograms of 
compensating fuel and a rotation factor of 0.9. This 
situation is compared to 100 kilograms of compensating fuel 
and a rotation factor of 1.0 in Figure 23. Only very slight 
differences occur between Figures 23 and 25. Mission life 
is a little longer in the retrograde orbits using a rotation 
factor of 0.9 and a little shorter in the prograde orbits 
using a rotation factor of 0.9. 

When the satellite orbit is lowered to an altitude at 
perigee of 200 kilometers (AO = 1.03) with 1000 kg of fuel, 
mission life is reduced in the retrograde-equatorial orbit 
below that of the prograde-equatorial orbit. The observed 
peaks at 70 and 140 degrees also shift leftward. The 
observed dimple seen in Figure 22 at the critical 
inclination of 63 degrees may be geopotential-induced or an 
artifact of the Kozai approximations as discussed by Taff 
[Ref. 2:pp 332-342] . 

A further explanation of the plot of mission life versus 
inclination shown in Figure 22 is presented here. A more 
detailed study of this input variable is .worthy of a 
complete paper itself. 

The low mission life found at low values of inclination 
is due to the combination of the point of perigee position 
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and oblate-atmosphere effects. The low inclination orbits 
travel through the "fat part” of the oblate atmosphere near 
the equator. As the inclination increases, the point of 
perigee moves away from the center of the atmospheric 
density-bulge and the satellite's orbit moves up out of the 
thick equatorial atmosphere. These two factors combine to 
produce the observed rapid rise in mission life seen in 
Figure 22. Beyond near-polar values of the inclination, the 
predicted mission life again falls. The mission life valley 
seen at 120 degrees is probably due to maximum orbital time 
near the atmospheric density-bulge center. 

Figure 23 shows the effect of starting with 100 
kilograms of initial drag-compensating fuel. Note the many 
peaks and valleys. To confirm that these peaks and valleys 
represent actual predicted mission life trends and not 
simply data scatter, the granularity was decreased to 1/100. 
Figure 24 is a plot of these results, and confirms the 
trends outlined in Figure 23. The many peaks and valleys 
superimposed over the general trend of mission life seen in 
Figure 22 are probably due to geopotential-induced orbital 
drift. 

The results of the sensitivity and trend analysis of the 
orbital inclination XI have been discussed with reference to 
in Figures 22 through 25. Future investigation of the 
plotted data may reveal more insight toward understanding 
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the underlying laechanisms . The sensitivity analysis plots 
show no trends that appear unreasonable. 

The sensitivity analysis for the input variables 
effecting the orbital-perigee atmospheric density-bulge 
encounter geometry (category one) have now been presented. 
The next set of analyses are those of the second category, 
which contains the input variables affecting the slowing 
power of the atmosphere. These variables are: decimetric 
solar flux F107, drag coefficient CD, cross-sectional area 
A, specific impulse SPI, and the atmospheric rotation factor 
D. The initial mass of drag-compensating fuel is also 
included in this category because it improves satellite 
mission-life resistance to atmospheric drag. 

The real-world expected effect on mission life as a 
result of varying the decimetric solar flux F107 is to vary 
the heating effectiveness of the sun. As the decimetric 
solar flux increases, the density of the atmosphere will 
increase . 

The Jacchia J60 model atmosphere used in the propellant 
longevity model subroutine JAC60 uses the F107 flux as a 
multiplier for the exponential density term. The expected 
model predictions of mission life as the F107 flux is varied 
through its range of 50 to 3 00 is a multiplicatively 
decreasing plot. 

The results of the sensitivity analysis for predicted 
mission life as the decimetric solar flux is varied are 
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presented in Figures 26 and 27. Figure 26 is the result of 
10 kg of initial maneuvering fuel, and Figure 27 is the 
result with 100 kg of fuel. The predicted mission life in 
Figure 2 6 ranges from 24.99 to 3.96 days. The 100 kg fuel 
case in Figure 27 predicts a mission life ranging from 
270.12 to 43.07 days. These results are in agreement with 
expectations . 

The next set of analyses consider effects of the drag 
coefficient, satellite cross-sectional area, and specific 
impulse. All of these factors are multipliers acting to 
change the resistance of the satellite to atmospheric drag. 
The drag coefficient CD is a proportionality constant in the 
drag force model. The satellite cross-sectional area A 
changes the size of the satellite that the atmosphere 
actually acts against. As satellite cross-sectional area 
increases it will take multiplicatively more energy to 
compensate for the orbital energy lost by the satellite 
working against the drag force. The specific impulse SPI is 
a multiplier that expresses the effectiveness of the 
expelled fuel mass in overcoming the drag force. 

The rotation factor D is the ratio of the atmospheric 
rotation rate to the earth's rotation rate. Its effect is 
to change the satellite to atmospheric relative velocity as 
the inclination changes. A logical range of D is between 
the values 0.9 and 1.0. The sensitivity analysis range of 
0.1 to 10.0 was chosen to show model sensitivity for a wide 
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Figure 26. Mission Life versus F10.7 
Flux (10 kg of Fuel) 
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Figure 27. Mission Life versus F10.7 
Flux (100 kg of Fuel) 
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range of D. As presented in Chapter II Equation 12, the D 
factor is combined with the inclination and the ratio of the 
radius-at-perigee to the velocity-at -perigee producing the 
6 factor. This factor is used to convert the satellite 
tangential velocity into the satellite-to-atmospheric 
relative velocity. 

The propellant longevity model uses A, CD, 6 , and SPI 
as proportionality constants, multiplicatively changing the 
effectiveness of the expelled fuel mass. The routine PLEP 
combines these four factors into a single constant DEL given 
by: 



DEL = A * CD * 6 * 3.14159/SPI (14) 

This constant DEL is passed to the routine KOZAI which in 
turn passes it to the routine DRAG. In the routine DRAG, 
DEL is used as the final multiplier in determining required 
fuel for each orbit. 

The expected model -predicted mission-life as a result of 
changing these factors is a multiplicative increase when A, 
CD, or 6 are increased linearly. A multiplicative decrease 
should occur when SPI is increased linearly. 

Expected results of varying D from 0.1 to 10.0 in the 
near-polar retrograde orbit are shown in Figure 28. Note 
that a very slight lowering of mission life is seen when the 
rotation factor increases. Figures 29 and 30 show the 
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Figure 28. Mission Life versus Rotation Parameter 
(10 kg of Fuel and Default Inclination) 
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Figure 29. Mission Life versus Rotation Parameter 
(10 kg of Fuel and Prograde Orbit) 
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Figure 30. Mission .Life versus Rotation Parameter 
(10 kg of Fuel and Retrograde Orbit) 
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plotted result of varying D through the same range of 0.1 to 
10.0. However, in Figure 29 an inclination of 25 degrees is 
used instead of the default 95 degrees, and in Figure 30 an 
inclination of 155 degrees is used. In Figure 29 it is seen 
that the mission life multiplicatively increases as D is 
linearly increased. Figure 30 shows a multiplicative 
decrease in mission life as D is linearly increased. These 
trends are due to the change in the satellite-to-atmosphere 
relative velocity. The increasing mission life expectancy 
within a strongly prograde orbit (Figure 29) is a result of 
the decreasing satellite-to-atmosphere relative velocity. 
The decreasing mission life expectancy within a strongly 
retrograde orbit (Figure 30) is a result of increasing 
satellite-to-atmosphere relative velocity. The plots shown 
in Figures 28, 29 and 30 are in agreement with expected 
results. 

The results of the sensitivity analysis plots on mission 
life as a result of varying CD, A, and SPI are shown in 
Figures 31 through 35. 

Figures 31 and 32 result from varying the satellite 
drag coefficient CD from 1.0 to 3.5. Figure 31 represents 
the 10 kg fuel case and Figure 32 the 100 kg case. The 
predicted mission life range in Figure 31 is 31.18 to 8.63 
days, and in Figure 32 from 308 to 99.58 days. 

The sensitivity analysis mission-life trend plots of 
satellite cross-sectional area A are shown in Figures 33 and 
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Figure 31. Mission Life versus Drag 

Coefficient (10 kg of Fuel) 
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Figure 32. Mission Life versus Drag 

Coefficient (100 kg of Fuel) 
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Figure 33. Mission Life versus Cross-Sectional 
Area (10 kg of Fuel) 



110 



MISSION LIFE (DAYS) 

100 200 300 400 500 600 700 




Figure 34. Mission Life versus Cross-Sectional 
Area (100 kg of Fuel) 
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Figure 35. Mission Life versus Specific 
Impulse (10 kg of Fuel) 
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34. Figure 33 is the 10 kg fuel case and Figure 34 the 100 
kg case. In both cases, A was varied from 0.000010 to 
0.0001000 square kilometers. The ranges mapped by these 
plots are 80 to 0.64 days and 736.29 to 108.7 days, 
respectively . 

Figure 35 is the analysis of motor specific impulse SPI 
varying from 200 to 3000 with 10 kg of initial fuel. The 
apparent linear shape of the SPI plot is due to its inverse 
multiplicative action and the range of the constant DEL (14) 
over which it acts. As SPI increases, the drag compensating 
fuel is more effective. Linearly increasing SPI has the 
same effect as starting at the default value * on the cross- 
sectional plots and moving along the plot to the left. Note 
that the shape of cross-sectional mission life plot in this 
range is closer to the expected mirror image of the SPI 
plot. 

As noted earlier, mission life plots showing "knees” 
have significance to research and design studies. The knees 
in the cross-sectional area mission life plots are very 
pronounced. Less distinctive knees are seen in the mission- 
life plots of the decimetric solar flux and the drag 
coefficient. The default value of .the cross-sectional area 
is right in the knee range of the predicted mission-life 
plots, as seen in Figures 34 and 35. 

The final input variable in the second category is the 
initial propellant mass W. The expected mission life of a 
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drag compensating satellite as fuel is increased linearly 
should be a linear increase in the predicted mission life. 
A given amount of fuel produces a certain mission life. 
Twice that amount of fuel in the same general configuration 
should produce about twice the mission life. The expected 
shape of the sensitivity analysis plot of predicted mission 
life for linearly increasing initial propellant mass is a 
approximately straight line. Figure 36 is the result of 
varying W on mission life. The plot is approximates a 
straight line as expected. The ripples are probably due to 
geopotential-induced orbital-drift or to movement of the 
atmospheric density-bulge during mission life. 

The results of the sensitivity and trend analysis 
mission-life plots for CD, A, D, SPI, and W agree with 
expectations. All input variables in the second category 
have been examined. The next input variables to be examined 
are in the third category. 

Consider now the input variables that affecting the 
depth of penetration of the satellite into the atmosphere at 
perigee. The variables that affect the satellite altitude- 
at perigee are the semi-major axis and the orbital 
eccentricity. These variables show the most pronounced 
effects on mission life. The real-world effect on mission 
life when linearly changing the altitude at perigee is an 
exponential increase in mission life. The effect of varying 
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Figure 36. Mission Life versus 

Initial Propellant Mass 
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the eccentricity and the semi-major axis on the altitude at 
perigee is summarized by 

rp = a * (1 - e) (15) 

Here, rp is the radius at perigee, a the semi-major axis, 
and e the eccentricity. A linear increase in the semi- 
major axis results in a modified exponential increase in 
mission life. A linear increase of the eccentricity results 
in a modified exponential decrease in mission life. 

The propellant longevity model computes the satellite 
altitude at perigee in the subroutine SALT. The subroutine 
DRAG then passes the altitude at perigee to the routine 
JAC60 where it is used to develop an exponent for an 
empirical model of atmospheric density. 

Mission life results of varying the orbital eccentricity 
ES from 0.01 to 0.1 are presented in Figures 37 and 38. 
Figure 37 is the 10 kg fuel case, and Figure 38 the 100 kg 
fuel case. The predicted mission life ranges from 14 to 0 
days in Figure 37, and 150 to 0 days in Figure 38. 

An analysis of the semi-major axis AO is presented in 
Figure 39 with a range of 1.025 to 1.175 earth radii. In 
this case, 1 kg of fuel is allocated initially. The 
resulting mission life is from 0 to 384 days. 

The plots of mission life resulting from varying the 
semi-major axis AO and the orbital eccentricity ES, show 
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Figure 37. Mission Life versus Orbit 

Eccentricity (10 kg of Fuel) 
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Figure 38. Mission Life versus Orbit 

Eccentricity (100 kg of Fuel) 
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Figure 39. 



Mission Life versus Semi-Major 
Axis (1 kg of Fuel) 
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major knees in their curves. The * symbol indicates the 
default value of the varied element. Note that the amount 
of fuel used in the analysis of the semi-major axis is 1 kg. 
The use of the usual 10 kg produces a plot that exceeding 
the reasonable range of mission life. The results shown in 
Figures 37 through 39 are in agreement with expectations. 

The input variables in the fourth category (no effect on 
mission life) are the mean anomaly ETU, initial satellite 
weight WT, average decimetric solar flux FEAR, and the 
geomagnetic index AKP. FEAR and AKP have no effect on 
mission life due to the limitations of the Jacchia J60 
atmosphere model. In the real world these two factors may 
combine with the decimetric solar flux F107 to produce 
changes in the atmospheric density of three orders of 
magnitude. 

The model theory used by Dr. Parks [Ref. 1] to predict 
the propellant longevity of a satellite doing intrack micro- 
thrusting to overcome atmospheric drag is independent of the 
satellite mass. The fuel required to overcome drag-induced 
orbital energy loss is based strictly on the satellite's 
shape and size. 

Figure 40 presents the sensitivity analysis plot of 
predicted mission life as a result of varying the satellite 
initial mass WT. The plot is a straight horizontal line, as 
expected . 
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Figure 40. Mission Life versus Satellite 
Mass (10 kg of Fuel) 
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The final analysis is concerned with the mean anomaly. 
In the near-circular default orbit, the mean anomaly and the 
eccentric anomaly are nearly the same. The Keplerian mean 
anomaly functions to position the satellite in its 
designated orbit. For a mission life lasting several days, 
it should make no difference where the satellite is 
initially positioned in its orbit. Rather, it's the orbit 
itself that makes the difference. 

Analysis plots of the m.ean anomaly ETU are presented in 
Figures 41 and 42. The incremented range is from 0 to 179 
degrees. Figure 41 is the 10 kg fuel case and Figure 42 is 
the 100 kg fuel case. Predicted results range from 13.74 to 
13.87 days in Figure 41, and from 149 to 154 days in Figure 
42. 

The results in Figure 41 are as expected, an approximate 
horizontal line. In Figure 42, a valley is seen in the 40 
to 120 degree range. The valley is only four days deep, but 
is difficult to explain. Further research may reveal its 
underlying mechanism. 

The results of the trend and sensitivity analysis have 
been presented. Mission life has been found to be most 
sensitive to those factors that affecting the orbital radius 
at perigee. The main unexpected results are found in the 
analysis of the inclination and the mean anomaly. The 
pattern of many peaks and valleys found in the 100 kilogram 
fuel case (Figures 23, 24 and 25) warrants further research. 
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Figure 42. Mission Life versus Mean 
Anomaly (100 kg of Fuel) 
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The four day dimple seen in the mean anomaly plot (Figure 

1 ^. 

41) also warrants further study. 





V. SUMMARY 



This final chapter presents the summary and conclusions 
of the thesis. Recommendations for model use are also 
given. 

The low-earth-orbit satellite propellant longevity model 
has been tested in a sensitivity analysis. Based on the 
results of the sensitivity analysis the model reasonableness 
is confirmed. The first chapter presented a model overview, 
model motivation, and the process employed to verify the 
model. Chapter I also discusses some changes made to the 
model computer program for the IBM 3033 at the Naval 
Postgraduate School. Chapter II presents the orbital 
mechanics and atmospheric concepts as background for model 
understanding and it's operation. In Chapter III the model 
and its computer program were discussed. Operating 
instructions specific for the computer program use at the 
Naval Postgraduate School are also presented there. The 
sensitivity analysis results are presented in Chapter IV. 

The propellant longevity model predicts how long the 
fuel of a low-earth-orbital satellite thrusting to overcome 
atmospheric drag will last before starting uncontrolled 
decay into the atmosphere. Model verification reasonable- 
ness tests reveal no inconsistencies with the established 
theory. The results of the sensitivity analysis are 
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therefore in agreement with the expectations. The 

irregularities in the model predictions (occurring with 
inclinations of less than 22.5 degrees) are due to computer 
round off error, not the model itself. 

The propellant longevity model currently uses the 
Jacchia J60 model atmosphere. This atmosphere model may 

significantly reduce the real-world accuracy of predicted 
mission lifetimes. 

The short term confidence factor in predicted mission 
life might be increased by the employment of the Jacchia J77 
model atmosphere. One drawback to the use of the Jacchia 
J77 model atmosphere is that it may double the computer run 
time. 

One solution to the trade-off of increased accuracy of 
the Jacchia J77 model versus the shorter computation time of 
the Jacchia J60 model, would be to allow the user to select 
an atmosphere option. 

Transportation of the computer program of the model to 
an IBM AT could provide more flexibility in preliminary 
design analysis. 

This model is a tool for the low-earth-orbit satellite 
system planner in designing and managing satellite systems 
that have precise orbital element control requirements. 
These systems are those with SDI and reconnaissance 
missions . 
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This model, with the addition of the Jacchia J77 model 
atmosphere and a subroutine to calculate fuel required to 
compensate for aspherical geopotential forces, could be used 
to optimize SDI defensive satellite constellation design. 

Existing real-world data to compare with the model 
predictions is limited to a very classified system. The 
final step in verifying the propellant longevity model would 
be to make this comparison test. 
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APPENDIX A 



COMPUTER PROGRAM LISTING 



xxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxx 



XX XX 

XX PROPELLENT LONGEVITY 

FOR 

LOW-ALTITUDE EARTH ORBIT SATELLITES 
(A.D. PARKS) 

- 






5^ THIS PROGRAM IS THE CODING OF THE THEORY CONTAINED IN NSWC X 

^ TR 83-243. IT WAS USED TO GENERATE THE NUMERICAL EXAMPLES ON 5C 

PAGES 13 TO 17. INITIALLY CODED FOR THE CDC 6700 IN FORTRAN 5^ 
^ IV BY A.D. PARKS, THE CODE HAS BEEN TRANSPORTED TO THE IBM 3033 
yi AND TRANSLATED TO FORTRAN 77 BY LT. C.D. NOBLE. ADDITIONAL DOC- ^ 

X UMENTATION WAS ADDED DURING THE TRANSPORTATION AND TRANSLATION. 5^ 

yi yi 

X THE REFERENCE MATERIAL CITED IN THE CODE IS MOSTLY FORM THE x 

)€ NAVAL SURFACE WEAPONS CENTER (NSWC) TECHNICAL REPORTS (TR). ^ 

X NSWC TR # ARE CITED IN THE CODE AS COMMENTS. OTHER REFERENCES X 

yi ARE SITED BY NUMBER. THEY ARE: ^ 

yi yi 

yi REFERENCE 1 BROUWER, D., "SOLUTION OF THE PROBLEM OF yi 

yi ARTIFICIAL SATELLITE THEORY WITHOUT DRAG", yi 

yi THE ASTRONOMICAL JOURNAL, VOL 64, NO 1274, X 

X PP. 378-397, 1959. ^ 

yi yi 

yi REFERENCE 2 AGRAWAL, B.N., "DESIGN OF GEOSYNCHRONOUS 5^ 

^ SPACECRAFT", PRENTICE-HALL, 1986 . 

yi yi 

yi REFERENCE 3 JACCHIA, L.G., "THERMOSPHERIC TEMPERATURE, ^ 

yi DENSITY AND COMPOSITION: NEW MODELS”, x 

yi SPECIAL REPORT 375, SMITHSONIAN ASTRO- X 

X PHYSICAL OBSERVATORY, 1977. x 

yi yi 

X REFERENCE 4 JACCHIA, L.G., "EMPIRICAL MODELS OF THE x 

X THERMOSPHERE AND REQUIREMENTS FOR IMP- X 

X ROVEMENTS", AFGL-TR-81-0195, ADA101946, x 

X 1981, 

yi yi 

X REFERENCE 5 JACCHIA, L.G., "ANALYSIS OF DATA FOR THE X 

X DEVELOPMENT OF DENSITY AND COMPOSITION x 

X MODELS OF THE UPPER ATMOSPHERE”, AFGL-TR- x 

X 81-0230, AD106420, 1981. X 

X yi 

X REFERENCE 6 ORR, L.H., "ORBITAL LIFETIME PROGRAM", x 

X NASA TM 87587, 1985. x 

X yi 

yi REFERENCE 7 TAFF, L.G., "CELESTIAL MECHANICS: A COMP- x 

X UTATIONAL GUIDE FOR THE PRACTIOMER", X 

X JOHN WILEY & SONS. INC., 1985. x 

X yi 

X REFERENCE 8 FLEAGLE, R.G., "AN INTRODUCTION TO ATMOS- X 

X PHERIC PHYSICS", ACADEMIC PRESS, 1980. X 

X ^ 
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MODULE 



PLEP 



SNALYS 



KOZAI 



BRAUER 



POSVEL 



NWTRPH 



FUNCTIONS 

XSPA-XSPR 



MEAN 



GEOP 



FUNCTION 



SYSTEM INITIALIZATION AND CONTROL. 

CALLED BY: CALLS: 

(SUBROUTINES) TIMES CALLED 
SNALYS 1 

KOZAI 1/ITERATION 

SENSITIVITY MODE SCANNED ELEMENT GENERATOR. 
CALLED BY: CALLS: 

(SUBROUTINES) TIMES CALLED 
PLEP xxxx xxxx 

SYSTEM WORKHORSE, DRIVER AND RECORD KEEPER. 



NSWC TR 83-135. 
CALLED BY: 

PLEP 



CALLS: 

(SUBROUTINES) 

BRAUER 

PERIOD 

MEAN 

THRST 

GEOP 

DRAG 



TIMES CALLED 
1 
1 
1 

1/(SAT.REV. ) 
1/(SAT.REV. ) 
1/(SAT.REV. ) 



CONVERTS NAVSPASUR BROUWER MEAN ORBITAL ELEMENT 
SET INTO A BROUWER OSCULATING ORBITAL SET. 

A CODING OF THE BROUWER-LYDDANE THEORY. 

NSWC TR 83-107, 83-135, BROUWER (REF. 1). 



CALLED BY: 
KOZAI 



CALLS: 

(SUBROUTINES) 

xxxx 



TIMES CALLED 
xxxx 



SOLOC 



COMPUTES SATELLITES POSITION AND VELOCITY IN 
INERTIAL SPACE FROM THE OSCULATING ORBITAL 
ELEMENTS. NSWC TR 82-387 PP. 2A-25. 

CALLED BY: CALLS: 

(SUBROUTINES) TIMES CALLED 
SATLOC NWTRPH 1 

A NEWTON-RAPHSON SOUUTION TO KEPLERS EQUATION. 
NSWC TR 82 387 P 23. 

CALLED BY: CALLS: 

(SUBROUTINES) TIMES CALLED 
POSVEL xxxx xxxx 

COMPUTES SHORT PERIODIC EFFECTS ON ORBITAL 
ELEMENTS. CALLED BY SUBROUTINE MEAN. NSWC 81-A56 
PAGES 9-11 AND NSWC TR 83 293 REF 3. 

CALLED BY: CALLS: 

(SUBROUTINES) TIMES CALLED 
MEAN xxxx xxxx 

USES THE WALTER METHOD TO CONVERT OSCULATING 
ORBITAL ELEMENTS TO KOZAI LIKE MEAN ELEMENTS 
USING FIRST ORDER PERIODIC VARIATIONS. 

SIMILAR TO THE METHOD IN NSWC TR 83 135 PP 16-18. 
CALLED BY: CALLS: 

(FUNCTIONS) TIMES CALLED 
KOZAI XSPA-XSPR 1 

COMPUTES AVERAGED SECULAR CHANGES IN KOZAI MEAN 
ELEMENTS INDUCED BY GRAVITATIONAL FIELD THROUGH 
FOURTH ZONAL HARMONIC (J2-J9 TERMS). NSWC TR 81- 
956. EQS.(2.5). 

CALLED BY: CALLS: 

(SUBROUTINES) TIMES CALLED 
KOZAI xxxx xxxx 

COMPUTES THE RIGHT ASCENSION AND THE DECLINATION 
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X 




OF THE SUN FOR 


THE MODEL ATMOSPHERE 


. NSWC TR 81- 


X 


X 




456 PP. 20-21. 






X 


X 




CALLED BY: 


CALLS: 




X 








(SUBROUTINES) 


TIMES CALLED 


X 


¥ 




DRAG 




xxxx 


X 


X 


PERIOD 


COMPUTES THE ANOMALISTIC PERIOD. 




X 

X 


3C 




CALLED BY: 


CALLS: 




X 


X 






(SUBROUTINES) 


TIMES CALLED 


X 


X 




KOZAI 


xxxx 


xxxx 


X 



X DRAG 

X 

X 

X 

X 

X 

X 

X 

X 

X 

X 

X 

X 



X 



COMPUTES DRAG 


EFFECTS AND FUEL MASS 


DECREMENT 


X 


REQUIRED TO MAINTAIN THE SEMIMAJOR 


AXIS VIA 


X 


IN TRACK MICROTHRUSTING. NSWC TR 82 


-243 EQ. 


X 


58. THIS ONE 


EQUATION IS THE HEART 


OF THE MODEL. 


X 


CALLED BY: 


CALLS: 




X 




(SUBROUTINES) 


TIMES CALLED 


X 


KOZAI 


SOLOC 


1 


X 




SATLOC 


3 


X 




SALT 


1 


X 




JAC60 


3 


X 




(FUNCTIONS) 




X 




MMBSIR 


1 


X 




XSPM 


1 


X 



X 



X JAC60 

X 

X 

X 

X 

X 

X 



COMPUTES ATMOSPHERIC DENSITY BASED ON THE JACCHIA x 
(60) MODEL ATMOSPHERE. CALLED THREE TIMES TO x 
COMPUTE DIFFERENT VALUES FOR THE DENSITY FOR EACH ^ 
REVOLUTION: MAX, MIN AND OP. X 
CALLED BY: CALLS: x 

(SUBROUTINES) TIMES CALLED ^ 
DRAG xxxx xxxx x 



X 

X THRST 

X 

X 

X 

X 

X 

X 

X 

X SALT 

X 

X 

X 

X 

X 



X 

CALCULATES CHANGES TO KOZAI ORBITAL ELEMENTS IF x 
AN ORBITAL ADJUST IS SCHEDULED. CALCULATES THE x 
MASS OF FUEL USED TO PERFORM THE ADJUST. x 

NSHC TR 83-31 EQ . 16, AND NSWC TR 81-456 EQ . 2-26 x 

CALLED BY: CALLS: x 

(SUBROUTINES) TIMES CALLED x 

KOZAI xxxx xxxx X 

X 

COMPUTES SATELLITE ALTITUDE. NSWC TR 81-456 EQ . X 
4.34. X 

CALLED BY: CALLS: x 

(SUBROUTINES) TIMES CALLED x 

DRAG xxxx xxxx x 

X 



X SATLOC 

X 

X 

X 

X 



COMPUTES SATELLITE LOCATION. RIGHT ASCENSION, x 
DECLINATION AND ITS GEOMAGNETIC LATITUDE. x 
CALLED BY: CALLS: x 

(SUBROUTINES) TIMES CALLED x 
DRAG POSVEL 1 x 



X 



X 



X FUNCTION 

X (MMBSIR) 

X 

X 

X 

X 

X 



COMPUTES BESSEL FUNCTIONS. AN IMSL ROUTINE x 

RESIDENT ON THE IBM 3800. THIS WILL VARY WITH X 

LOCAL INSTALLATIONS. x 

CALLED BY: CALLS: X 

(SUBROUTINES) TIMES CALLED X 
DRAG xxxx xxxx x 

X 



xxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxx 



X 

X 

X 

X 

X 

X 

X 

X 

X 

X 



xxxxxxxxxxxxxxxxxxxxxxxxxxxxxx 

xxxxxxxxxxxxxxxxxxxxxxxxxxxxxx 

XX XX 

XX PROGRAM PLEP XX 

XX 5(X 

xxxxxxxxxxxxxxxxxxxxxxxxxxxxxx 

xxxxxxxxxxxxxxxxxxxxxxxxxxxxxx 
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PROGRAM PLEP 



THIS MODULE IS THE INPUT, CONTROL AND INITIALIZATION SECTION FOR THEx 
5< REST OF THE SYSTEM. IT ACCEPTS THE *FIVE CARD' FORMAT OF NAVSPASUR X 
X AS THE BROUNER MEAN ORIBITAL ELEMENT SET. THE SPECIFIC SATELLITE X 
3^ PARAMETERS ARE ENTERED. ATMOSPHERIC MODEL ELEMENTS ARE INITIALIZED. ^ 
3^ FINALLY, SYSTEM CONTROL CONSTANTS ARE SET TO GENERATE THE OUTPUT. A 3^ 
3^ FURTHER DESCRIPTION OF THE ABOVE FUNCTIONS IS CONTAINED IN THE ^ 
^ DESCRIPTION OF ARGUMENTS. 3t 

yi 



yk X 

3^ DESCRIPTION OF ARGUMENTS 3^ 

yi X 

> INPUT ARGUMENTS 3^ 

yi yi 

yi PROGRAM CONTROLS 3^ 

yi yi 

3^ ITERl INTEGER, MAXIMUM NUMBER OF REVOLUTIONS DESIRED 3t 

y^ FOR CONSIDERATION. A CONTROL PARAMETER. 3( 

3^ ITER2 INTEGER, ORBITAL PRINT FREQUENCY. A CONTROL 3( 

3^ PARAMATER. SETTING ITER2 = TO 100 WILL REDUCE 3^ 

3^ TABULAR OUTPUT BY ONLY PRINTING OUTPUT EVERY 100 3^ 

3^ ORBITS 3^ 

X PSENS INTEGER, CONTROL PARAMETER. PSENS=1 SPECIFIES * 

X SENSITIVITY ANALYAIS MODE. PSENS=0 SPECIFIES x 

X STANDARD PROPELLENT LONGIVITY MODE. X 

X lET INTEGER, CONTROL PARAMETER. IET=1 PRINTS THE X 



X INITIAL ANOMALISTIC PERIOD. IET=0 SUPRESSES PRINT. x 

X- SYSTEM RESET TO 0 WHEN SENSITIVITY MODE IS SELECTED 



X NOA INTEGER, NUMBER OF SCHEDULED ORBIT ADJUSTS. THIS X 

X IS A CONTROL PARAMETER WITH ALLOWED VALUES FROM x 

X ZERO TO TEN. THE ONLY ORBIT ADJUSTS CURRENTLY x 

X SUPPORTED ARE CHANGES TO THE SEMI-MAJOR AXIS. x 

X IRV(K) INTEGER ARRAY. (K=10 MAX). THE ELEMENTS OF THE X 

X ARRAY SPECIFY THE REVOLUTION NUMBER DURING WHICH X 

X EACH OF THE ORBITAL ADJUSTS ARE PERFORMED. K MUST x 

X EQUAL NOA. X 

X DA(K) INTEGER ARRAY. (K=10 MAX) THE ELEMENTS OF THE X 

X ARRAY SPECIFY THE CHANGE IN KILOMETERS OF THE X 

X SEMIMAJOR AXIS. K MUST EQUAL NOA. X 

X PSIZE INTEGER, CONTROL PARAMETER THE NUMBER OF INTERVALS X 

X IN THE SENSITIVITY ANALYSIS. X 

X X 



X SATELLITE PHYSICAL PARAMETERS X 

^ X 

X SPI REAL, MOTOR SPECIFIC IMPULSE IN SECONDS. AGRAWAL x 

X (REF. 2:PP. 163-166) GIVES SAMPLE VALUES. X 

X WT REAL, INITIAL MASS OF THE SATELLITE IN KILOGRAMS, x 

X (INCLUDING THE FUEL MASS) x 

X W REAL, INITIAL MASS OF FUEL IN KILOGRAMS. x 

X CD REAL, DRAG COEFFICIENT OF THE SATELLITE. x 

X DIMENSIONLESS. NOR.'IALLY VARIES FROM ONE TO x 

X AROUND THREE, DEPENDING ON VELOCITY, DENSITY, x 

X COMPOSITION OF ATMOSPHERE AND PHYSICAL DESIGN. x 

X IS MOST OFTEN SET AT 2.20 FOR ALTITUDES BETWEEN X 

X 200 AND AOO KM. FOR ALTITUDES OF 1000 KM CD IS X 

X SET AT JUST ABOVE 3.0. JACCHIA (REF. 3:P 2). x 

X A REAL, SATELLITE CROSSESTIONAL AREA IN SQUARE x 

X KILOMETERS. X 

X X 



X ATMOSPHERIC SPECIFIC PARAMETERS X 

X X 

X D REAL, RATIO OF ATMOSPHERIC TO EARTH ANGULAR x 

X ROTATION RATES. A CONSTANT IN EQUATION 11 OF X 



X NSWC TR 83-243. UNITLESS. THE ACTUAL VALUE'WILL X 

X DEPEND ON LATTITUDE AND ALTITUDE. RANGE .8 TO 1.0. X 

X F107 REAL, DECIMETRIC SOLAR FLUX. JACCHIA (REF 485). x 

X THIS VALUE IS USED BY THE INDUSTRY TO REPRESENT x 

X THE XUV FLUX IN THE THERMOSPHERE. THE XUV FLUX x 
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X IS THE FACTOR CREATING THE SUN FACING DIURNAL x 

X DENSITY AND TEMPERATURE BULGE IN THE THERMOSPHERE . X 

X IT IS DIRECTLY RELATED TO SOLAR ACTIVITY. IT X 

X VARIES WITH THE 27 DAY SLOAR ROTATION AND 11 YEAR X 

X SUNSPOT ACTIVITY CYCLE. ORR (REF. 6;APPENDIX B) X 

X INDICATES RANGES OF 73.3 TO 2A2.9. XUV FLUX AND X 

X THE GEOMAGNETIC HEATING EFFECTS CAN CAUSE X 

X DENSITY TO CHANGE BY THREE ORDERS OF MAGNITUDE AT x 

X 600 KILOMETERS, THE 10.7 CM FLUX I“S DECTABLE FROM X 

X THE GROUND AND THEORETICALLY PARELLELS THE XUV x 

X FLUX. THE XUV FLUX IS NOT MEASUREABLE FROM THE X 

X GROUND DUE TO NEAR COMPLETE UPPER ATMOSPHERE X 

X ABSORPTION. X 

X FBAR REAL, THE AVERAGE 10.7 CM FLUX. THE EFFECTS OF THEx 

X XUV HEATING DO NOT RETURN TO THE ORIGINAL UNHEATEDX 

X CONDITION IN AN EARTH ROTATION. A HIGH HEATING X 

X CONDITION WILL HAVE SOME LONGER TERM EFFECTS. x 

X- THESE LONGER TERM EFFECTS ARE MODELED HERE BY x 

TAKING A MOVING 90 DAY 10.7 CM AVERAGE. x 

X FBAR RANGES FROM A LOW OF 73 TO A HIGH OF 230. x 

X AKP REAL, GEOMAGNETIC ACTIVITY INDEX, JACCHIA (REF. A X 

X AND 5). GEOMAGNETIC ACTIVITY ALSO CAUSES HEATING X 

X IN THE THERMOSPHERE. THE EXACT MECHANISM IS NOT X 

X UNDERSTOOD. IDEALLY QUEIT MAGNETIC CONDITIONS x 

X CORRESPOND TO AN INDEX OF 0. MAXIMUM ACTIVITY FOR X 

X WHICH SIGNIFICANT DATA EXISTS ARE RELATED TO AN X 

X INDEX OF BETWEEN 6 AND 7 . NORMAL ACTIVITY x 

X INDICATES AN INDEX OF BETWEEN l.AND 2. X 

X TVE REAL, TIME OF VERNAL EQUINOX PASSAGE. USED TO* x 

X TRANSFORM TIME FROM INERTIAL TO EARTH FIXED x 

X REFERENCE FRAMES. FOUND IN THE ASTRONOMICAL ALMANAC 

X FOR THE YEAR OF CHOICE. (MODIFIED JULIAN DAYS). x 

X A DESCRIPTION IS IN TAFF (REF. 7:PP. 103-10A). x 

X X 

X NAVSPASUR ORBITAL ELEMENTS x 

X X 

X THE NAVSPASUR ORBITAL ELEMENT SET IS THE BROUWER MEAN X 

X ELEMENT SET. BROUWER (REF. 1) IS THE BEST DESCRIPTION FOR x 

X THIS ELEMENT SET. x 

X X 

X UJD REAL, TIME OF EPOCH IN MODIFIED JULIAN DAYS. X 

X THIS IS THE TO IN THE EQUATION FOR THE MEAN ANOMALY 

X M. ( M = MO + N X (T - TO) ). TAFF (REF. 7). x 

X ETU REAL, INITIAL MEAN ANOMALY MO. MEASURED IN DEGREESX 

X HAS A RANGE OF 0 TO 360. X 

X GO REAL, MEAN ARGUMENT OF PERIGEE. RANGES FROM 0 TO x 

X 360 DEGREES. NOT DEFINED FOR CIRCULAR ORBITS. x 

X HO REAL, MEAN RIGHT ASCENSION OF THE ASCENDING NODE. X 

X RANGES FROM 0 TO 360 DEGREES. NOT DEFINED FOR X 

X EQUITORIAL ORBITS X 

X ES REAL, MEAN ECCENTRIY. FOR THIS THEORY RANGES FROM X 

X 0.01 TO 0.1. THESE ARE LOW ECCENTRICITY ORBITS. x 



X 

X 

X 

X 

X 

X 

X 

X 

X 

X 

X 

X 

X 

X 

X 

X 

X 

X 

X 



XI REAL, MEAN INCLINATION. RANGES FROM 0 TO 180 X 

DEGREES. VALUES OF 0 AND 130 ARE IN EQUITORIAL x 
ORBITS. INCLINATION OF 90 IS A POLAR ORBIT. RANGESX 
OF 0 TO 90 ARE IN PROGRADE ORBITS. INCLINATIONS x 
90 TO 180 ARE IN RETROGRADE ORBITS. FOR AN X 

INCLINATION OF 0 OR 180 HO IS NOT DEFINED. THERE x 
IS NO ASCENDING NODE IN EQUITORIAL ORBITS. x 

RND REAL, RATE OF CHANGE OF MEAN MOTION. MEASURED IN x 

REVOLUTIONS PER (DAY SQUARED) . VALUES OF .00001 x 
AND ABOVE ARE CONSIDERED SIGNIFICANT WHEN X 

ATMOSPHERE INDUCED. REFERENCE (TAFF) P137 X 

ESD REAL, ECCENTRICITY DECAY RATE MEASURED IN x 

INVERSE SECONDS. X 

AO REAL, KAULA SEMIMAJOR AXIS. MEASURED IN x 

MEAN EQUITORIAL EARTH RADII (6378.165 KM). RANGES x 
FROM 1.0259A TO 1.16847 TO MAINTAIN AN ALTITUDE OFx 
100 TO 1000 KM FOR RADIUS AT PERIGEE. THIS MUST BEX 
CAREFULLY SCALED TO ECCENTRICITY TO PREVENT X 

CRASHING THE SATELLITE. PERIGEE ALTITUDES ABOVE X 
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600 KM ARE ESCENTIALLY AVOVE THE ATMOSPHERE. X 

^ ADOT REAL, SEMIMAJOR AXIS DECAY RATE. X 



X 

DIMENSION SNSARC20), SNSMINC20), IRV(IO), DACIO), SNSDEU20) 

X 

3( THE SET OADJ IS COMMON 'TO PLEP, KOZAI, THRST AND SNALYS. 

COMMON / OADJ/ NOA, IRV, DA, SPI, WT, PSENS 
X 

^ THE SET XXXX IS COMMON TO PLEP AND PERIOD. x 

X 

COMMON / XXXX / lET 
X 

5(9(5(5(X9<?^X9^M5^X9(XXXX5(X9(5<5<5(5(5(X5(9«9(5(X9f9(X9(X5<5(5(9(9«9(X9(XXXX5(X9(5(3(9(X9(XX3(9(9(X5(XX5^XXX5( 

3( THE SET SOLAR IS COMMON TO PLEP AND DRAG. 5( 

X 

COMMON / SOLAR / F107 , FBAR , AlCP , TVE 
INTEGER PSIZE, PSENS, PSYNS 
REAL MINSN 

X 

^ INITIALIZING CONSTANTS. ^ 

X 

DATA R2 / 541.15E-06 / ,DEG /57 . 295779513D0 / , ROE /637S.165D0/ 
X INPUTTING CONTROL PARAMETERS, 5( 

5(5(5(X5<3(X3(B(X9t5(XX3(3<3(XX9(X5(5(3(5(9C3(3(X9(5(5(3(B(9(3(XXXB(XXXXX3(3(XXXX5()(XXX5(^XX9(X5(9(XX5(XX3( 

READ X, ITERl, ITER2, NOA, PSENS, lET, PSIZE 

X 

^ INPUTTING SATELITE SPECIFIC PARAMETERS. 

X 

READ 5(, CD, A, D, SPI, W, WT 
X 

xxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxx 
X INPUTTING EARTH AND ATMOSPHERIC PARAMETERS. 

XXXXxXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXX 

X 

READ X , F107 , FBAR , AKP , TVE 

X 

xxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxx 
X INITIALIZING BROUWER MEAN ORBITAL ELEMENTS FROM THE NAVSPASUR * 
3^ FIVE CARD FORMAT ORBITAL DATA SET. ^ 

xxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxx 
X 

READ 5 , UJD , ETU , HO , GO , ES , XI 
5 FORMAT ( /8X,F14.8,5(1X,F8. A)) 

READ 10 , RND , ESD 
10 FORMAT ( 20X, F11.9, 19X,E12.5 ) 

READ 15 , AO , ADOT 
15 FORMAT ( /8X,F11.5,1X,E14.7) 

X 

xxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxx 
X INPUTTING THE REVOLUTION NUMBER AND CORRESPONDING CHANGE IN 
X SEMIMAJOR AXIS x 

xxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxx 
X 

IF ( NOA .EQ. 0) GO TO 30 
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DO 20 K = 1, NOA 

READ IRVCK) , DACK) 

20 CONTINUE 

5( THIS NEXT SERIES OF CODE SETS THE INITIAL VALUES FOR THE SCANNING 
^ ARRAY MATRIX SO AS TO ALLON THE RECOVERY OF THE ORIGINAL INPUT 5( 
^ ELEMENTS IF SENSITIVITY MODE IS NOT SELECTED. 

30 DO 90 ICN = 1,20 

SNSDELCICN) = 1.0 
SNSMINCICN) = 0.0 
SNSARCICrO = 0.0 
90 CONTINUE 

^ THE NEXT SERIES OF WRITE STATEMENTS OUTPUTS THE ORIGINAL INPUT 
•X PARAMETERS. UNITS ARE SUPPLIED FOR CLARIFICATION. x 

WRITE(6,98) 

WRITE(6,99) 

WRITE(6,103) 

WRITE(6,99) 

WRITE(6,104) D, F107, FBAR, AKP, TVE 
WRITE(6,10S) CD, A, SPI, W, WT 
WRITEC6,106) UJD, ETU, HO, GO, XI 
WRITE(6,107) RND, ES, ESD, AO, ADOT 
WRITE(6,103) ITERl, NOA 
IFCNOA .GT. 0) THEN 
WRITEC6,109) 

DO 27 NN = 1, NOA 

WRITE(6,110) IRV(NN), DA(NN) 

27 CONTINUE 

END IF 
WRITE(6,99) 

X THIS NEXT CODE SETS THE INTIAL INPUT PARAMETERS TO SOME RETAINERS 
X VARIABLES THAT ALLOW RETURNING TO THE ORIGINAL VALUES IF THE X 

X SENSITIVITY MODE IS SELECTED. THIS PREVENTS HAVING TO REREAD ALL x 
X THE INPUT PARAMETERS BACK IN FOR EACH RUN OF THE SENSITIVITY SCAN.X 

XXXXXX5(XXXXXXXX5(5(XX5(XXXXXXXXX5^X9^^XX5(5(XXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXX 

V1=CD 

V2=A 

V3 = D 

V4=SPI 

V5=W 

V6=WT 

V7=F107 

V3=FBAR 

V9=AKP 

V10=TVE 

V11=UJD 

V12=cTU 

V13=H0 

VIA^GO 

V15=ES 

V16=XI 

V17=RND 

V13=ESD 

V19=A0 

V20=AD0T 

XX5^X^XXXX5(XX5CXXXXXXXXX?(5CXXXXXXXXXXXXXXXXXXX^XXXXXXXXXXX5(XX5^9CXXXXXXXXXX 

X THIS IS THE ITERATIVE SCANNING ROUTINE. ARGUMENTS ARE INCREMENTED x 
X FRACTIONALLY THROUGH THEIR SPECIFIED RANGE AS SPECIFIED IN THE ^ 
X INPUT TO THE SUBROUTINE SNALYS (PSIZE ITERATIONS OF PLEP). X 

X IF SENSITIVITY MODE IS NOT SELECTED PSIZE = 1. 

x^xxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxx 
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IFCPSENS .EQ. 0) THEN 
PSIZE=0 
END IF 

DO 40 PSYNS = 0. PSI2E 

X THIS NEXT SET OF CODE IS ENVOKED WHEN PSENS IS EQUAL TO ONE. X 

X IT SETS UP FOR AND CALLS THE SENSITIVITY ARGUMENT ARRAY SNALYS FOR^ 
SUBSEQUENT SCANNING. CPSENS = 1 SELECTS SENSITIVITY MODE x 

IF (PSENS .EQ. 1) THEN 
IET = 0 

X THE SUBROUTINE SNALYS INTIALIZES THE SCANNED INPUT PARAMETERS TO x 
X BE INCLUDED IN THE STUDY. IT ALSO PRINTS THE SENSITIVITY MATRIX. 

IFCPSYNS .EQ. 0) THEN 

CALL SNALYSC PSIZE, SNSAR, SNSMIN, SNSDEL, SUMCK, DELSN, 
1 MINSN) 

ELSE 

CD = VI 

RECOVERY OF THE ORIGINAL INPUT PARAMETERS AFTER THE FIRST RUN OF 5^ 

9( A SENSITIVITY ANALYSIS. X 

X 

A = V2 
D = V3 
SPI= V4 
W = V5 
WT = V6 
F107= V7 
FBAR= V8 
AKP= V9 
TVE=V10 
UJD=V11 
ETU=V12 
HO =V13 
GO =V14 
ES =V15 
XI =V16 
RND=V17 
ESD=V18 
AO =V19 
ADOT=V20 
END IF 
END IF 

ASYNS=REAL(PSYNS) 

ASIZE=REAL(PSIZ5) 

IFCPSENS .EQ. 1) THEN 

IFCSUMCK .EQ. 1,0) THEN 

SCNFRC = ASYNS^DELSN + MINSN 

ELSE 

SCNFRC = ASYNS/ASIZE 
END IF 
END IF 

X THIS NEXT CODE CREATES THE CURRENT VALUE FOR THE SCANNED INPUTS. x 
X IF SENSITIVITY IS SELECTED THE PARAMETER WILL SCAN FROM MIN TO MAXX 
X WITH A GRANULARITY OF l/PSIZE PER SCAN 

CD = SNSARC1))^(SNSMIN(1)+ASYNSXSNSDEL(1))-CDX(SMSAR(1)-1 .0) 

A = SNSAR(2)X(SNSMIN(2)+ASYNSXSNSDEL(2) )-A x(SNSAR(2)-l .0) 

D = SNSAR(3)9^(SNSMIN(5)+ASYNSXSNSDELC3))-D x( SNSARC 3) -1 . 0 ) 

SPI= SNSAR(4)9^(SNSMIN(4)+ASYNS5^SNSDEL(4))-SPI5^CSNSAR(4)-1 .0) 
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W = SNSAR(5)X(SNSMIN(5)+ASYNSXSNSDEL(5) )-W SNSARC 5 ) -1 . 0 ) 

WT = SNSAR(6)^(SNSMIN(6)+ASYNS^SNSDEL(6) )-WT^(SNSAR(6)-l .0) 

FI 07= SNSARC7)3((SNSMIN(7)+ASYNS5(SNSDEL(7))-F107 5((SNSAR(7)-1 .0) 
FBAR= SNSAR(8)X(SNSMINC8) + ASYNS5(SNSDEL(8))-FBAR^((SNSAR(8)-1 .0) 

AKP= SNSAR(9)5<(SNSMIN(9)+ASYNSXSNSDEL(9))-AKP5<(SNSAR(9)-1.0) 

TVE = SNSAR(10)5((SNSMIN(10) + ASYNSXSNSDEL(10) )-TVE3((SNSAR(10)-l .0) 

UJD = SNSARC11)?((SNSMIN(11)+ASYNS5(SNSDEL(11))-UJDX(SNSAR(11)-1 .0) 
ETU=SNSAR(12)X(SNSMINC12)+ASYMSB(SNSDEL(12) )-ETU^( SNSARC 12)-1 . 0 ) 

HO =SNSAR(13)X(SNSMIN(13)+ASYNS^SNSD£L(13))-H0X(SNSAR(13)-1 .0) 

GO =SMSAR(14)X(SMSMINC14)+ASYNSXSNSDEL(14) )-G03<(SNSAR(14)-l ,0) 

ES =SNSAR(15)X(SNSMIN(15)+ASYNS5€SNSDELC15))-ES?((SNSAR(15)-1 .0) 

XI =SMSAR(16)X(SNSMIN(16)+ASYMS5<SNSDEL(16))-XIX(SNSAR(16)-1 .0) 
RND=SNSAR(17)X(SNSMIN(17)+ASYNSXSMSDEL(17))-RNDX(SMSAR(17)-1 .0) 

ESD = SNSAR(18)X(SNSMINC18) + ASYNS?(SNSDEL(18))-ESD3((SNSAR(18)-1 .0) 

AO =SNSAR(19)X(SNSMIN(19)+ASYNS^SNSDEL(19) )-A0X(SNSAR(19)-l .0) 
ADOT=SNSAR(20)^(SNSMIN(20)+ASYNS3(SNSDEL(20))-ADOTX(SNSAR(20)-1 .0) 

^ DEL IS THE CONSTANT MULTIPLIER FOR THE INTEGRAL IN EQUATION 33 OF ^ 

X NSWC TR 83-243 ^ 

DEL = CDXAX3.14159D0/SPI . 

X 

xxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxx 
^ IF THE SENSITIVITY ANALYSIS MODE IS NOT SELECTED, THEN THE SYSTEM ^ 

X PRINTS STANDARD MISSION LIFE TIME TABULAR DATA. ^ 

xxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxx 

X 

IF (PSENS .EQ. 0) THEN 
BCOF = CD ^ A 
WRITE (6,98) 

WRITE (6,99) 

WRITE (6,100) 

WRITE (6,99) 

WRITE (6,101) BCOF, SPI, W 
WRITE (6,102) F107, FBAR, TVE 

END IF 
X 

xxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxx 
^ CONVERTING DEGREES TO RADIANS 

xxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxx 

X 

XI = XI / DEG 
GO = GO / DEG 
HO = HO / DEG 
FLO" = ETU / DEG 

X 

xxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxx 
THIS NEXT CODE SEQUENCE CONDITIONS THE ELEMENTS AS PER NSWC TR 82-387 
^ PAGE 7. IT CONVERTS THE KAULA SEMIMAJOR AXIS AO IN EARTH RADII TO x 
^ BROUWER MEAN SEMIMAJOR AXIS' VIA EQNS 2.1 AND 2.2. THE SEMIMAJOR AXISX 
^ DECAY RATE ADOT IS CONVERTED IN A LIKE MANNER VIA 2.3, 2.4, AND 2.5 x 
xxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxx 

X 

FN = 1. / ( 1. - ES X ES )XX1.5D0 
TA = SIN ( XI ) 

TB = TA X TA 

TE = 1. - (3.5(TB) / 2. 

TAD = (( 3.)(R2) / (2.?^A0>^A0))xTE^FN 
AK = AO 

AO = (A05(((1. + 2.^TAD) / (1. - TAD) )xxo . 66666667D0 ) ^ROE 
XDOT = TAD)((3.^ESD.^(ES/(1.-ES^ES)) - 2 . ^( ADOT/AK) ) 

ADT = (ADOT/AK)XAO + 2 . 5eAK)eR0E^XD0T?(( ( 1 . + 2.^TAD)^ 

1 ((1. - TAD)XX 5))5(X(-0,3333333333D0) 

XXXXXXXXXXXXXXXXXXXXXXXXXXXKXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXX 

X TIME IS NOW CONVERTED FROM SECONDS TO DAYS ^ 

xxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxx 

ADT = ADT / 86400. DO 
ESD = ESD / 86400. DO 
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RND = RNDx2.D0)f3.14159D0X(l./(86400.D0))C)C2) 

X THE SUBROUTINE KOZAI IS THE OPERATING WORKHORSE FOR THE PROPELLENT x 
X LONGEVITY COMPUTATIONS. IT ALSO DECREMENTS PROPELLENT AS IT IS USEDX 

X 

CALL KOZAICADT.AO. ES,XI.FL0.G0,H0.UJD,W.DEL,ITER1.ITER2, 

1 SCNFRC, D) 

X 

xxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxx 
X FORMAT STATEMENTS USED IN THE PLEP INITIALIZATION PROGRAM. X 

xxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxx 

X 



98 FORMAT ( IHl ) 

99 F0RMATC//5X. ' •,//) 

100 F0RMATC//15X, 'PROPELLANT LIFE ESTIMATOR',//) 

101 F0RMATC//5X, »CDA=',E12.5, ' SQ KM '/5X, 'SPECIFIC IMPULSE^', 

XF7.2,' SEC/5X, 'INITIAL WEIGHT OF PROPELLANT^ ' , 

XF7.2,' KG') 

102 FORMAT (//5X, ' F107= ' , F8 . 2/5X, ' FBAR= ' , F8 . 2/5X, 

X ' TIME OF VERNAL EQUI NOX= ' , E22 . 1 A ) 

103 F0RMATC15X, 'INITIAL INPUT CONDITIONS') 

lOA FORMATC///1X, 'ATMOSPHERIC INPUT PARAMETERS ' //5X, ' D (RELATIVE ATMOS 
XPHERIC ROTATION) = ' , FI 0 . 4/5X, ' FI 07 (DECIMETRIC SOLAR FLUX) = ', 
XF7.2/5X, 'FBAR (90 DAY AVG. F107) = ' , F7 . 2/5X, ' AKP (GEOMAGNETIC IND 
XEX) = ',F6.3/5X, 'TVE (TIME OF VERNAL EQUINOX PASSAGE) = ',F16.9, 

X' MODIFIED JULIAN DAYS'///) 

105 FORMATdX, 'SATELLITE PHYSICAL PARAMETERS V/5X, ' CD (DRAG COEFFICIEN 
XT) = ',F5.2,' UNITLESS'/5X, 'A (CROSSECTI ONAL AREA) = ',F11.9,' SQU 
XARE KILOMETERS'/5X, 'SPI (MOTOR SPECIFIC IMPULSE) = ',F8.2,' SECOND 
XS'/5X,'W (INITIAL FUEL MASS) =« ',F9.2,' KI LOGRAMS '/5X, ' WT (INITIAL 
X SATELLITE MASS) = ',F11.2,' KILOGRAMS'///) 

106 FORMAT( IX, ' NAVSPASUR ORBITAL EL EMENTS ' //5X, ' UJ D (EPOCH) = ',F1A.8 
X,' MEAN JULIAN DAYS '/5X, ' ETU (MEAN ANOMALY) = ',F8.4,' DEGREES'/5X 
X,'H0 (MEAN RIGHT ASCENSION OF THE ASCENDING NODE) = ',F8.4,' DEGRE 
XES'/5X,'G0 (MEAN ARGUMENT OF PERIGEE) = ',F8.4,' DEGREES'/5X, 

X'XI (MEAN INCLINATION) = ',F8.4,' DEGREES') 

107 FORMAT(5X, 'RND (RATE OF CHANGE OF MEAN MOTION) = ',F11.9,' REVOLUT 
XTIONS PER DAY SQUARED'/ 

X5X,'ES (MEAN ECCENTRICITY) = ',F8.4,' UNITLESS '/5X, 'ESD (ECCENTRIC 
XITY DECAY RATE) = ',F16.9,' 1/SECONDS '/5X, 'AO (KAULA SEMI-MAJOR AX 
XIS) = ',F11.5,' MEAN EARTH RADII '/5X, ' ADOT (SEMI-MAJOR AXIS DECAY 
XRATE) = ',F16.9,' '///) 

108 FORMATdX, 'CONTROL INPUT PARAMETERS '//5X, 'MAXIMUM REVOLUTIONS PER 

, XITERATION ', 17 ,/ 5X ,' NUMBER OF SCHEDULED ORBIT ADJUSTS TO SEMI-MAJO 
XR AXIS ',12/) 

109 F0RMAT(//5X, 'SEMI-MAJOR AXIS INCREASE SCHEDULE'//5X, ' REV ADJUST', 
X5X, 'CHANGE IN KILOMETERS'/) 

110 F0RMAT(5X, I7,7X,F6 .2) 



X 



X 

X 

X 

X 

X 

X 

X 

X 

X 

X 

X 

X 

X 



40 CONTINUE 

:'0P 

END 



xxxxxxxxxxxxxxxxxxxxxxxxxxxxxx 

xxxxxxxxxxxxxxxxxxxxxxxxxxxxxx 

XX XX 

XX SUBROUTINE SNALYS 3(X 

XX . XX 

xxxxxxxxxxxxxxxxxxxxxxxxxxxxxx 

xxxxxxxxxxxxxxxxxxxxxxxxxxxxxx 



SUBROUTINE SNALYS (PSIZE, SNSAR, SNSMIN, SNSDEL, SUMCK, DELSN, 
IMINSN) 
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DIMENSION SNSAR(20), SNSMIN(20), SNSMAXC20), SNSDEU20) 

CHARACTER5(4 SNSCAR 
INTEGER N, PSIZE 

REAL SNSAR, SNSMIN, SNSMAX, SNSDEL, SUMCK, DELSN, MINSN 
WRITE (6,9S) 

WRITE(6,120) 

WRITE(6,100) 

WRITE(6,120) 

WRITE(6,200) PSIZE 
ASIZE=REAL(PSIZE) 

READ(5,102) 

SUMCIC = 0.0 
DO 20 N = 1, 20 

READ 101, SNSCAR, SNSAR(N), SNSMIN(N), SNSMAX(N) 

SNSDEL(N) = (SNSMAX(N)-SNSMIN(N))/ASIZE 
IF (SNSAR(N) .EQ. 1.0) THEN 
SUMCK = SUMCK +1.0 
DELSN = SNSDEL(N) 

MINSN = SNSMIN(N) 

WRITE(6,300) SNSCAR, SNSMIN(N), SNSMAX(N), 

1 SNSDEL(N) 

END IF 
20 CONTINUE 

IF (SUMCK .GT. 1.0) THEN 
DELSN = 0.0 
MINSN = 0.0 
END IF 

WRITE(6,120) 

WRITEC6, 400) 

98 FORMAT ( IHl ) 

120 F0RMATC//5X, ' •,//) 

100 F0RMATC//17X, 'SENSITIVITY ANALYSIS FOR ' , / 16X, ' L OW EARTH ORBIT •, 

X 'SATELLITES', /19X, 'PROPELLENT LONGEVITY',//) 

200 FORMAT (5X, ' TABLE OF INPUT PARAMETERS SCANNED FROM MINIMUM TO', 

X /5X,' MAXIMUM WITH THE INDICATED INCREMENT PER RUN. ',/5X, 

X ' THE GRANULARITY IS 1/ ' , 14, • . '/5X, ' ALL OTHER INPUT PARAMETERS 
X ARE RETURNED TO THEIR', /5X,' INITIAL VALUES AT EACH INCREMENT.' 

X ///'INPUT PARAMETER ',8X,' MINIMUM ',13X, 'MAXIMUM ', 13X, 'INCREMENT'/) 

102 FORMAT (IX) 

101 FORMAT ( A4, IX, F3.1, F16.9, 4X, F16.9) 

300 FORMAT (5X, A4, 6X, 3(3X, F16.9)) 

400 FORMAT(//, lOX, ' SENSITIVITY SCAN RESULTS FOR ABOVE SPECIFICS'// 
15X,'REV NUMBER', 8X, 'TIME(DAYS) ',15X, 'FRACTION OF SCAN') 

RETURN 

END 



A AA AA 

X SUBROUTINE KOZAI 

A AA AA 
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yK ^K ^K ^K ^ ^K ^K 



XXXXX3(XXXXXX3(XB(3(XXXX9(3(3(XXX3(3(3(X 



SUBROUTINE KOZAI IS THE SYSTEM WORKHOUSE. IT CALLS BRAUER TO CON-3( 
VERT THE BROUWER MEAN ORBITAL ELEMENTS OF THE NAVSPASUR INPUT TO X 
BROUWER OSCULATING ELEMENTS. THESE OSCULATING ELEMENTS ARE CON- x 
VERTED TO KOZAI-LIKE MEAN ELEMENTS IN THE SUBROUTINE MEAN. THE X 
SUBROUTINE PERIOD COMPUTES THE INITIAL ANOMALISTIC PERIOD. THE x 
SUBROUTINE GEOP COMPUTES THE EFFECTS OF THE OBLATE GEOID ON THE x 
ORBITAL ELEMENTS (THROUGH THE TERM). SUBROUTINE THRST COMPUTES x 
THE EFFECTS ON THE ELEMENTS DUE TO PERIGEE BURNS TO CHANGE X 
THE SEMIMAJOR AXIS DURING IRV(K) BY THE AMOUNT DACK), FINALLY, x 
SUBROUTINE DRAG COMPUTES THE THE AMOUNT OF FUEL NEEDED TO OVERCOME x 
THE DRAG EFFECTS BY CONTINUOUS INTRACK- MICROTHRUSTING . X 
KOZAI KEEPS TRACK OF THE FUEL USED AND ORBITAL ELEMENT CHANGES. X 
WHEN THE FUEL IS EXHAUSTED KOZAI RETURNS SYSTEM CONTROL TO PLEP. x 



SUBROUTINE KOZAI ( ADT, AO , ES, XI , FL 0 , GO , HO , UJD, W, DEL , ITERl , ITER2, 

1 SCNFRC D) 

DIMENSION X0SCC6) , XMC6), DXDTC6) , XDDC6), IRV(IO), DA(IO) 
DIMENSION DXMC6) 

COMMON / OADJ/ NOA, IRV, DA, SPI, WT, PSENS 

DATA B0,B2,B3,B4,B5,B6 / . 398603254E+06 , - . 1755528999E+11, 

1 .26386647738E+12, . 1 063073996E+16 , . 8056 05022E+18 , 

1 .7292115856E-0‘^/ 

DATA DEG/57 .29577951 3D0/,CN2/0 .0/,Al, El, RNl/0. 0,0. 0,0 .0/,DT/0 .0/ 
DATA PI / 3.1^159D0 / 

DATA NPLOT / 0 / 

IRC = 0 ‘ 

PCOUNT=0 
KK = 0 

PI2 = 2.D0X PI 

CALL BRAUER ( BO , B2, B3 , B4, B5 , DT, AO , ES, XI , FLO , GO , HO , CN2, A, CE, Cl , CL , 
1 G,H,A1,E1,RN1 ) 



IF (PSENS. EQ. 0) THEN 

WRITE(6,100) ITERl 

100 F0RMAT(//5X, 'MAXIMUM NUMBER OF ORBITAL REVOLUTIONS= ' , II 0 ) 
WRITE(6,101) ITER2 

101 F0RMAT(//5X, 'DATA PRINTED EVERY ',13,' ORBITAL REVOLUTIONS') 
SCNFRC =0.0 

END IF 

CALL PERIOD ( AO, ES, XI, PA ) 

XOSC(l) = A 
X0SC(2) = CE 
X0SC(3) = Cl 
XOSC(^) = G 
X03C(5) = H 
XQSC(6) = CL 
CALL MEAN (XOSC,XM) 

ILINE = 0 
CIW = XM(3)xDEG 
GW = XM(^)XDEG 
HW = XM(5)XDEG 
CLW = XM(6)XDEG 
IF (PSENS .EQ. 0) THEN 
WRITE(6,102) 

WRITE(6,103) (XM(K),K=1,2),CIW,GW,HW,CLW 
END IF 

103 F0RMAT(/5X, 'SEMIMAJOR AXIS= ' , E22 . 1<^ , ' KM'/ 

X 5X, 'ECCENTRICITY=',E22.14/ 

X 5X, 'INCLINATICN=',E22.1^, ' DEG'/ 

X 5X, 'ARGUMENT OF PERIGEE= ' , E22 . 1^ , ' DEG'/ 

X 5X, 'ASCENDING NODE= ' , E22 . 1<^, ' DEG'/ 

X 5X,'MEAN AN0MALY=',E22.1^, ' DEG') 
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102 F0RMATC//5X, *K0ZAI MEAN ELEMENT SET’) 

DELT = PA 

CALL THRST ( XM, DXM, IRC, KK, PA, DWOA) 

DO 10 1=1, ITERl 

T = ( I - 1 ) K DELT / 86A00. 

TT = UJD + T 

15 CALL GEOP (XM,DXDT) 

16 DO 20 JJ = 1,6 

XM(JJ) = XM(JJ) + DXDT(JJ)B^DELT + DXMCJJ) 
IF ( JJ . LT . A ) GO TO 20 

XMC JJ ) = AMOD ( XM(JJ), PI2 ) 

20 CONTINUE 

IRC = 0 

DO 30 II = 1, NOA 

IF ( I.EQ. IRVCID) IRC = 1 
IF ( IRC .EQ. 1) KK = II 
30 CONTINUE 

CALL THRSTCXM, DXM, IRC, KK, PA, DMOA) 

X 

X CALCULATING THE DELI IN NSWC TR 83-2A3 EQN 11. 
RP=XM(1)J^(1.0-XM(2)) 

VP=SQRT(CB0/XM(1))X((1 .0+XM(2))/(l .0-XMC2)))) 
DEL1 = (1 .0-(RP/VP)5^D5^B6XCOS(XM(3)))3^3(2 
DDELT = DEL5(DEL1 



5^ 

X 

X 

X 

yi 



CALL DRAG ( XM, DDELT, TT, DELMG ) 

W = N + DELMG - DNOA 
WT = NT + DELMG - DWOA 
IF (PSENS .EQ. 0) THEN 

IF (PCOUNT.NE. (ITER2-D) THEN 
PCOUNT=PCOUNT+l 
GO TO 10 

ELSE 

PCOUNT=0 
END IF 

IF ( ILINE .GE. 50 ) ILINE = 0 
ILINE = ILINE + 1 

IF ( ILINE .EQ. 1 ) WRITE ( 6,10A) 

IF ( ILINE .EQ. 1) WRITE ( 6,105) 

WRITE ( 6, 106) I, TT, T, W 
END IF 

IF ( W.LE. 0.00) THEN 

IF (NPLOT .EQ. 1) THEN 
PRINT T, SCNFRC 

ELSE 

WRITE(6,107) I, T, SCNFRC 
END IF 
RETURN 
END IF 
10 CONTINUE 

IF (NPLOT .EQ. 1) THEN 
PRINT K, T, SCNFRC 

ELSE 

WRITE(6,107) I, T, SCNFRC 
END IF 

lOA FORMAT(lHl) 

105 F0RMATC/5X, ’REV NUMBER ’, 3X, ’ TIME(MJD) * ,8X, ’TIME(DAYS) ’, 8X, 
X ’PROPELLANT REMAININGC KG) ’ , // ) 

106 FORMAT(5X,I10,3X,F9.2,3X,F10.2,15X,F8.2) 

107 FORMAT(5X,I10, 8X, FI 0 . 2 , 15X, F16 . 9 ) 

RETURN 

END 



yiyiyiyiyiyiyiyiyiyiy^yiy^yiyiyiyiyiyiyiyiiiyiyiyiyiyiyiyiyi 

yiyiyiyiyiyiyiyiyi'Ayiyiyiyiyiyiy^yiyiyiyiyiyi'Ayiyf^y^yiyiyi 

yiyi 

SUBROUTINE BRAUER 
yiyi y^yi 
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SUBROUTINE BRAUERCBO, B2, B3, B4, B5, DT, A2P, E2P> CI2P, CL02P,G02P, H02P, 
1CN2,A,CE,CI,CL,G,H,ADT,RND,ESD) 
lOOA A2P2=A2P'^5$2 

A2P4 = A?P?5t3(? 

CN0=SQRT(B0/(A2P25(A2P)) 

E2P2 = E2P?^5(2 
ETA=3QRT(1 .-E2P2) 

SINEI=SINCCI2P) 

THETA=COS(CI2P) 

THETA2 = THETA'^X2 
THETA4 = TNETA2^5€2 
THETA6=THETA45^THETA2 
CJ2 = -B2/(2.?^BOXA2P2) 

ETA2=ETAXX2 
ETA3 = ETA25^ETA 
ETAA = ETA23(5^2 
CJ21P=CJ2/ ETA4 

CJ31P = 33/(B0^A2P25eA2P^ETAA5eETA2) 

CJ-ilP = (3.5(BA)/(8.)^B0XA2P^^ETA45^ETA4) 

CJ51P = B5/(B0)€A2P4XA2P^ETA4‘)^5<25(ETA2) 

FUN1 = 3.5^THETA2-1 . 

FUN2 = 1 .-5.5(TH£TA2 

SINEI2=SINEIXX2 

A1=A2PXCJ25(FUN1 

AO=-Al/ETA3 

A2 = 3.5(A2PXCJ2XSINEI2 

FUN5=1 .-11 .^THETA2-(40.XTHETA4)/FUN2 

FUN6= -FUN1-C8 .XTHETA4)/FUN2 

FUN4=THETA2/SINEI2 

FUN22 = FUN2X5^2 

r J21 P2 = r J71 P5^5^2 

E01P = ((E2P5^ETA2)X(3.5(CJ21P25(FUN5-I0.5e 
1CJ41P5€FUN6))/(24.5^CJ21P) 

E21P=-2.^E01P 

E31P=((35.?cCJ51PxE 2P2XETA25^SINEI)5€(FUN2-(16 .xTHETA4)/FUN2))/(96 .X 
1CJ21P) 

E11P = -.7 5D0)^E31P+(C .25D0)eETA2XSINEI)5e(CJ31P+.3125D0XCJ51P 
1?^(4.+3.XE2P2)X 

2( 1 . -9 . 5^THETA2-(24 . 5^THETA4 ) /FUN2) ) )/CJ21P 
CI0 = -(E2P5(THETA)/(ETA25(SINEI) 

CI2=CJ21P5<THETA5(SINEIK1 .5DO 
CI1 = E2P)^CI2X. 6666666 7 DO 
FUN7 = (-.5D05^ETA35^CJ21P)/E2P 

CL21P=(ETA3/CJ21P)5<( . 25D03^CJ21P25^FUN5- . 83333333D05(CJ41P5(FUN6 ) 
CL12P=CNOX(l .+1 .5D05(CJ21P5^£TA5(FUN1+.09375D05(CJ21P25^ETA)^(-15.+16 . 
lxETA+25 . XETA2+C30 . -96 . )^ETA-90 . XETA2 ) 5(THETA2+( 1 05 . +144 . xETA+25 . 
2XETA2)5€THETA4)+.9 375D05^CJ41PXETA5€E2P25€(3.-30 .xTHETA2+35.5€THETA4) ) 
CL22P=.5D0^CN05^CN2 

G21P = (1 ./(24.'^CJ21P))^(-3.5(CJ21P2X(2.+E2P2-11 .X(2.+3.B(E2P2)^THETA2 
1-40 . 2 . +5 . XE2P2) 5^THETA4/FUN2-40 0 . 5(E2P25^THETA6/FUN22 ) + l 0 . ?^CJ41?5^ 

2(2. 

3+£2P2-3.)^(2.+5.5CE2P2)xTHETA2-8 .5^(2.+5.55E2P2)5^THETA4>FUN2-8 0.)^E2P2X 
4THETA6/FUN22)) 

G12P=CN0X(-1 . 5D0XCJ21P5^FUN2+ . 09 375D05(CJ21P25€ ( -35 . +24 . 5^ETA+25 . 
1)^ETA2+ 

2(90 .-192. ^^ETA-126 .5eETA2) 5€THETA2+( 385. +36 0. 5(ETA+45.XETA2)5(THETA4) + 
3.3125D05^CJ41P5((21 .-9 .5^ETA2+(-27 0 . + 126 .XETA2)XTHETA2+(385.-189. 
4^ETA2)?$THETA4)) 

H2=l .5D05^CJ21P5(THETA 
H3 = -2.5€H2 

H1 = .666 66 667D05^E2P5(H2 

H31P = ((35.?<CJ51PXE2P25(E2P5^THETA)/(144.5^CJ21P))X( .5D0/SINED^(FUN2-( 
116 .5^THETA4)/FUN2)+SINEl5^(5. + (32.5(THETA2)/FUN2+80 .5^THETA4/FUN22) ) 
H11P=-.25X 

1 H31P+(( .25D05(E2P)^THETA)/(CJ21P5(SINEI))5((CJ31P+.3125D0^CJ51P5e 

2(4.+3.5(E2P2)5^(1 ,-9.5(THETA2-(24.XTHETA4)/FUN2) + l .875D0 5€CJ51P 
35(SINEI25^ 

4(4.+3.^E2P2)5((3. + (16 .XTHETA2)/FUN2+(40 .XTHETA4)/FUN22) ) 
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H21P = (E2P2)^THETA)/a2.^CJ21P)x(-3.)^CJ21P2)^(ll . + (30 .)^THETA2)/FUN2+ 
1(20 0.5^THETAA)/FUN22)+10 .)^CJ^1P)^(3. + (16 .3^THETA2)/FUN2+(A0 o^^THETAA)/ 
2FUN22)) 

H12P=CN05(TH£TAX(-3.XCJ21P+.375D0)^CJ21P2X(-5. + 12.)^ETA+9.5^ETA2+ 
lC-35.- 

236 .XETA-5.)CETA2))^THETA2) + 1 .25D0XCJ<4lPX(5.-3.3^ETA2)X(3.-7 .5CTHETA2)) 
AID=CJ51P/CJ21P 
AID2 = FUN2-(16 .)^THETA^)/FUN2 
C1 = 35./38A .)^AID)^ETA3)^E2P)^SINEI5^AID2 
AID3=THETA2/SINEI 
AIDA = THETA25<SINEI 
E2P3=E2P2^E2P 

C2 = 35./1152.)^AID)^((-E2P)^$INEIX(3.+2.XE2P2) + E2P33^AID3))^AID2+ 
12.5CE2P3)^AIDAX(5. + (32.)^THETA2)/FUN2+(30 .)^THETAA)/FUN22)) 

C3=l .-9.)eTHETA2-(2A.xTHETA<4)/FUN2 
AID5=CJ31P/CJ21P 

C^=.25D05(AID55€(-E2P)^AID3)+5./6A,XAIDX(-E2P)^AID35^(A. + 3.5CE2P2) + 
1E2P)^SINEI)^(26 .+9 .^<E2P2))5CC3-15./32.)^AIDxE2P5^AIDA)^(A.+3.)^E2P2)3^ 

2(3. + (16 .XTHETA2)/FUM2+(A0 .3^THETAA)/FUN22) 

C5 = E2P/(1 .+ETA3))^(3.-E2P23^(3.-E2P2)) 

C6 = (E2PX(-32.+31.x(E 2?2^E2P2)))/((A.+3.)^E2P23+ETA3C(A.+9 .3CE2P2) ) 
C7=.25D0)^AID5)^SINEl5^C5+5./64.3^C33eAID)^ETA2)^SINElxC6 
C8 = -.25D03<AID55CETA33^SINEI-5./64.)^AIDxETA3)^SINEI)e(A. + 9.3(E2P2))^C3 
0510 T=DT 

CL2P = CL12P^DT+CL22P3(DT)^9^2+CL02P + RND)^DT3(DT 
CL2P = AM0D(CL2P, 6 .233135307 17 96 DO) 

IF(CL2P)520,530,530 
520 CL2P=CL2P+6.2331353071796DO 
530 G2P = G12P)^DT+G02P 
H2P=H12PXDT+H02P 
SINEG=SIN(G2P) 

COSING=COS(G2P) 

D1E=SINEG5C(SINEG3^(E31P)^SINEG+E21P) + E11P) + E01P 
H1P = ((H31P5CSINEG+H21P))^SINEG+H11P)3^C0SING+H2P 

GPLP = G2P+CL2P+.5D03^(CL21P+G21P)xSIN(2 .XG2P) + (C1+C2)5CC0S(3 .3^ G2P) 

1 + (CA+C7)3^C0SING 
CL1P=CL2P 
U=CL2P 

100 DELTAU=(U-E2P3(SINCU)-CL2P)/(1 .-E2P3eCOS(U)) 

U=U-DELTAU 

I FCABSCDELTAU)-!. £-10)20 0,100 ,100 
20 0 U = U-(U-E2P3^SIN(U)-CL2P)/(1 .-E2P3^COS(U) ) 

E=U 

SINE1P=SIN(E) 

C0SE1P=C0S(E) 

G1P=G2P 

ADIVR=1./(1.-E2P3<C0SE1P) 

SINF1P=ADIVRXETAXSINE1P 
C0SF1P = ADIVR3^(C0SE1P-E2P) 

F1P=ATAN2(SINF1P,C0SF1P) 

I F(ABS(F1P-CL2P)-3.1A15926 53539300)220, 21 0,210 
210 STOP 

220 FUfJ3 = (l.+CN23^T)X3C.6666666 7D0 
C0SFG=C0S(2.X(G1?+F1P)) 

SINFG = SIN(2.3^(G1P+F1P)) 

ADIVR2 = ADIVR36^<2 
ADIVR3 = ADIVR>^>:3 

CI=CI2P+CI036D1E+CI13CSINF1P3^SINFG+(2 .3^CI13(C0SF1P+CI2)3^C0SFG 
FUfi3 = FlP-CLlP+E2P3€SINFlP 

3013 H = H1P+(2.3^H1hCO$F1P + H2)xSINFG-H13^SINF1P3^COSFG+H33^FUN3 
KFUN=H/6 .2331353071796D0 
FUN9=KFUN 

H=H-FUN9X6 .2831353071796D0 
IF(H)3022, 8023, 8023 
3022 H=H+6 .2831353071796D0 

8023 A = A2P/FUN3+AO + (A1+A2XCOSFG)3^ADIVR + ADT3^DT 
AID6=ADIVR2XETA2+ADIVR 
AID7-SIN(2.)^G2P+F1P) 

AID8 = SIN(2.3CG2P+3.3^F1P) 

D1=.25D0XCJ21P3^(6 . x( 5 . 3^THETA2-1 , )3^FUN3+( 3 . -5 . 3^THETA2)3^( 3 , 3^$INFG+ 
13.3(E2P3(AID7 +E2P3<AID8 )) 

D2=.25D03^CJ21P3^(2.3^(3.5^THETA2-1 , )3^(AID6 + 1 . ) 3^SINFlP+3 . x( 1 . -THETA2 )3€ 
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1 ( (-AID6+1 . )xAID7+(AID6+.33333333D0)xaID8)) 

AID9=C0S(2.xG2P+F1P) 

AID10=C0S(2.XG2P + 3.5^F1P) 

D3=-ETA2X.5DCXCJ21P^(1 .-THETA2)X(3.XAID9+AID10) 
ETA6I=1./(ETA3XETA3) 

D«^ = ETA6I?$(C5 +C0SF1PX(3. + E2P5<C0SF1PX(3.+E2P3^C0SF1P))) 

D5 = ETA6IX(E2P+C0SF1PX(3. + E2PXC0SF1P?^(3.+E2PXC0SF1P) )) 

D6= ETA23^CJ25(.5D0X( (3.XTHETA2-1 . ) )^DA+3 . 1 . -THETA2)XD5XC0SFG) + D3 

GAL=GPLP + D1 + (E2P5«ETA2)/(1.+ETA)XD2 

CE = (E2P-1.)X(1.+CN2XDT)XX.66666666666667D0 + 1.+D1E+D6 + ESD^^DT 
EDL=.5D05^E2P5(CL21P^SIN ( 2 . XG2P)+C8XCOSING+E2P«1XCOS (3.^G2P)- 
1 ETA3XD2 
AID1^=SIN(CL2P) 

AID15=C0S(CL2P) 

ESL=CE XAID14+EDL5(AID15 
ECL=CE XAID15-EDLXAID14 
CE=SQRT(ECL?^ECL + ESLXESL) 

CL=ATAN2(ESL,ECL) 

G=GAL-CL 

G=AM0D(G,6 .2831853071796D0) 

IF(G)3024, 3025, 8025 
3024 G = G + 6 .2831353071796D0- 
8025 RETURN 
END 
X 
X 

X xxxxxxxxxxxxxxxxxxxxxxxxxxxxxx 

X xxxxxxxxxxxxxxxxxxxxxxxxxxxxxx 

X XX XX 

X XX SUBROUTINE POSVEL xx 

X XX XX 

X xxxxxxxxxxxxxxxxxxxxxxxxxxxxxx 

X xxxxxxxxxxxxxxxxxxxxxxxxxxxxxx 

X 



X 



X 

X 

X 

X 

X 

X 

X 

X 

X 

X 

X 



SUBROUTINE POSVEL ( A , CE, Cl , CL , G, H, BO , R1 , R2, R3, VI , V2, V3, R, VEL ) 
CALL NWTRPH ( CL . CE/ E/ SINEP , COSEP) 

HSIN=SIN(H) 

HCOS=COS(H) 

GSIN=SIN(G) 

GCOS=COS(G) 

CISIN=SIN(CI) 

CICOS=COS(CI) 

A11=HCOSXGCOS-HSINXCICOSxgSIN 

A12 = -HCOSXGSIN-HSIN5^CICOS5^GCOS 

A21 = HSIN3(GCOS+HCOSXCICOSXGSIN 

A22=HC0SXCIC0SXGC0S-HSINXGSIN 

A31=CISIN5^GSIN 

A32=CISINXGC0S 

FUN=SQRT(1.-CEXCE) 

Rl = AX(Alix(COSEP-CE)+A125^(FUNxSINEP) ) 

R2 = AX(A215€(C0SEP-CE) + A22X(FUNXSINEP)) 

R3 = AX(A315^(C0SEP-CE)+A32X(FUNXSINEP) ) 

R=AX(1.-CE5«C0SEP) 

FUM1 = (SQRT(S05^A) )/R 

V1=FUN1k(A11X(-SIMEP)+A12XFUNxC0SEP) 

V2 = FUN1X(A21X(-SINEP)+A225(FUN?<C0SEP) 

V3 = FUtll^(A3lx(-SinE?)+A32XFUNXC0SEP) 

VEL=FUN1XSQRT(1 .-(CE^C0SEP)XX2) 

RETURN 

END 



xxxxxxxxxxxxxxxxxxxxxxxxxxxxxx 

xxxxxxxxxxxxxxxxxxxxxxxxxxxxxx 

XX XX 

5^x SUBROUTINE NWTRPH xx. 

XX XX 

xxxxxxxxxxxxxxxxxxxxxxxxxxxxxx 

xxxxxxxxxxxxxxxxxxxxxxxxxxxxxx 
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SUBROUTINE NWTRPH ( QL , B , E2, SE, CE) 

NEWTON RAPHSON SOLUTION TO KEPLERS EQUATION SUBROUTINE 
El=QL + (B5^SIN(QU)/(l.-B^COS(QL)) 

SE=SIN(E1) 

CE=COS(El) 

E2=£1 + (QL + B5(SE-E1)/(1 .-B^CE) 

IF(ABS(E2-ED- 1 . E-8 
E1 = E2 
GO TO 2 
SE=SINCE2) 

CE=C0S(E2) 

RETURN 

END 



XX XX 

XX FUNCTION XSPA xx 

XX XX 

xxxxxxxxxxxxxxxxxxxxxxxxxxxxxx 

xxxxxxxxxxxxxxxxxxxxxxxxxxxxxx 



THE FOLLOWING FUNCTIONS DETERMINE THE FIRST ORDER SHORT 
PERIODIC VARIATIONS AS A FUNCTION OF THE MEAN ELEMENTS. 

FUNCTION XSPA(XM) 

DIMENSION XMC6) 

REAL J(^) 

DATA J,R/1 .,0. 001 08276 DO, -0.00 00 0255D0, -0.000 00156 DO, 6378. 165D0/ 

F = XM(6)+(2.xXM(2)-XM(2)XX3./4+5./96 .XXM(2)XX5. )XSIN(XM(6)) 

1+(1 .25D0XXM(2)xx2.-11./24.xXM(2)xxA.+17 ,/192 . xxM( 2)xx 6 . )xSIN(2.x 
2XM( 

26)) + (13./12.xXM(2)xx2.-<^3./192.xXM(2)xx5. )XSIN(3.XXM(6)) 

P = XMCDXCl .-XM(2)XX2. ) 

RR = P/(l .+XM(2)XC0S(F)) 

XSPA= J(2)XRXR/XM(1)X((XM(1)/ RR )XX3 . X( ( 1 . -3 ./2 . XSINC XM( 3) )X$IN 
l(XM(3)))+3./2.xSIN(XM(3))xSIN(XM(3))xC0S(2.x(XMC4)+F)))- 
2(1 .-3./2.xsIN(XM(3))xSIN(XM(3)))x(l .-XM(2)XX2. )XX(-1 .5D0)) 

RETURN 

END 

X 

X 

X xxxxxxxxxxxxxxxxxxxxxxxxxxxxxx 

X xxxxxxxxxxxxxxxxxxxxxxxxxxxxxx 

X XX XX 

^ 5(X FUNCTION XSPE XX 

X XX XX 

X xxxxxxxxxxxxxxxxxxxxxxxxxxxxxx 

X xxxxxxxxxxxxxxxxxxxxxxxxxxxxxx 

X 

X 

FUNCTION XSPE(XM) 

DIMENSION XM(6) 

REAL J(A) 

DATA J,R/1. ,0.00108276D0,-0.00000255D0,-0.00000156D0,6378.165D0/ 
r = XM(6) + C2.^XM(2)-XM(2)*x3./4+5./96.xXN(2)5(X5. )5^SIN(XM(6) ) 

1 + (1.25D0xxM(2)b^x2.-11 ./24 . xxMC 2)xx4 . + 17 . / 1 92 . xXM( 2 ) xx6 . )XSIN(2.5^ 
2XM( 

36) ) + (13./12.xXM(2)Btx2.-43./192.xXM(2)5^x5. )*SIN(3.XXM(6)) 

P = XM(1)5C(1 .-XM(2)5fX2. ) 

XSPE = J(2)/2.x(R/P)XX2,X(1.-1.5D0xSIN(XM(3))xSIN(XM( 3))) 
lx((l .+1 .5D0XXM(2)XXM(2)-(1 .-XM(2)5^XM(2))5f^l,5D0)/XM(2) + 3.5((l .+ 
2XM(2) 

3xXfK2VA. )XCOS(F) + l .5DOXXM(2)3(COS(2.?^F) + XM(2)XXM(2)/4.*COS(3.XF)) 
4+3./8 .xJ(2)^(R/P)xx2.3(SIN(XM(3))xSIN(XM(3) )B (((1 .+11 ./A . XXM( 2)xx2 . ) 

5 ^C0S(2.^XM(A)+F)+XM(2)*XM(2)/4,XC0S(2.xXM(4)-F)+5.xXM(2)x 

6 C0S(2.5^XM(A) + FX2. ) + (7 . + 17 . /4 . XXM( 2) xXM( 2) )/3 . XCOSC 2 . 5CXM( A)+3 . *F) 

7 + 1 .5DOxXMC2)5^COS(2.5(XM(<:»)+4.xF)+XM(2)xXMC2)/4,xCOS(2.xXM(A)+5.xF) 

8 + 1.5DOxXM(2)^COS(2.5«XM(A)) ) 
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RETURN 

END 



X 

X 

3( 

X 

X 

X 

X 

X 



X 

X 

X 

X 

X 

X 

X 

X 

X 

X 

X 



xxxxxxxxxxxxxxxxxxxxxxxxxxxxxx 

xxxxxxxxxxxxxxxxxxxxxxxxxxxxxx 

XX XX 

FUNCTION XSPI XX 

XX XX 

xxxxxxxxxxxxxxxxxxxxxxxxxxxxxx 

xxxxxxxxxxxxxxxxxxxxxxxxxxxxxx 



FUNCTION XSPKXM) 

DIMENSION XM(6) 

REAL J(4) 

DATA J,R/1. ,0.00108276D0,-0.00000255D0,-0. 00000156D0,6378.165D0/ 
F = XM(6)+(2.xXM(2)-XM(2)xX3./4+5./96.xXM(2)xx5.)XSIN(XM(6)) 

1 + (1 .25D0XXM(2)XX2.-11./24.3(XM(2)XX4.+17 ./192 . XXM(2)XX6 . )XSIN(2.X 
2XM( 

36 )) + ( 13. /12.XXMC 2) XX2. -43. /192.XXMC 2 ) XX5. )XSIN( 3. XXM(6)) 

P = XM(1)X(1.-XM(2)XX2.) 

XSPI = 3./8.XJ(2)X(R/P)XX2.XSIN(2.5^XM(3))X(XM(2)xC0S(2.XXM(4) + F) 
1 +C0S(2.X(XM(4) + F))+XM(2V3.XC0S(2.XXM(4) + 3.XF)) 

RETURN 

END 



xxxxxxxxxxxxxxxxxxxxxxxxxxxxxx 

xxxxxxxxxxxxxxxxxxxxxxxxxxxxxx 

XX XX 

XX FUNCTION XSPW xx 

XX XX 

xxxxxxxxxxxxxxxxxxxxxxxxxxxxxx 
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FUNCTION XSPW(XM) 

DIMENSION XM(6) 

REAL J(4) 

DATA J,R/1 . ,0.001 08276D0,-0 .00 000255D0, -0.00000 156D0, 6378. 165D0/ 

F = XM(6)+(2.xXM(2)-XM(2)xx3./4+5./96.XXM(2)xx5.)XSIN(XM(6)) 

1 + (1.25D0XXM(2)3CX2.-11 ./24 . XXM( 2) XX4 . +17 ./192 . ^XM( 2) XX6 . ) XSINC 2 . 
2XM( 

26)) + (13./12.5(XM(2)XX2.-43./192.XXM(2))^X5. )XSIN(3.XXM(6)) 

P = XM(1)X(1.-XM(2)5<X2.) 

XSPW = 3./4.xJ(2)x(R/P)xx2.x(4.-5.xSIN(XM( 3))XSIN(XM(3))) 
1X(F-XM(6)+XM(2)XSIN(F)) + 1.5D0XJ(2)X(R/P)xx2.3((l ,-l .5D0XSIN(XM(3)) 
2XX2. ) 

3X((l,-XM(2)XXM(2)/4, )/XM(2)XSIN(F)+.5D0XSIN(2.xF)+XM(2)3(SIN(3.xF) 
4/12 

5.)-1.5D0xj(2)^(R/P)XX2.X( (SINCXMC 3) )XSIN( XMC 3) )/4 . +XM( 2) XX2 ./2 . X( 



61.-15./8 .XSIN(XM(3))XX2. ))/XM(2)xsiN(2.5(XM(4) + F)+XM(2)/16 .XSIN 
7(XM(3))':<5(2.XSIN(2.XXM(4)-F)+.5D0>^(1.-2.5D0XSIN(XM(3) )XX2.JXSIN(2.3( 

8(XM(4) + F) )-(7 ./12.5^SIN(XM(3))5^5<2.-XM(2)5<XM(2)/6 .X(l .-19 ./8 .)fSIN 
9(XM(3))^H2. ))/XMC2)3^SIN(2.>CXM(4) + 3.?<F)-3./3.3^SIN(XM(3))XX2.^SIN 
11(2.5SXM(4)+4.XF)-XM(2)/16 .3<SIN(XM(3) )XX2.5^SIN(2.XXM(4) + 5.3^F)) 
12-9./16 .XJ(2)5((R/P)XX2.3^SIN(XM(3))XX2.3^SIN(2.5^XM(4)) 

RETURN 

END 

X 
X 
X 
X 
X 
X 
X 
X 
X 
X 
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XX XX 

XX FUNCTION XSPO xx 

XX XX 

xxxxxxxxxxxxxxxxxxxxxxxxxxxxxx 

xxxxxxxxxxxxxxxxxxxxxxxxxxxxxx 
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FUNCTION XSPOCXM) 

DIMENSION XMC6) 

REAL J(4) 

DATA J,R/1., 0.00108276 DO, -0. 0000 0255D0, -0.00000156 DO, 6 378. 165D0/ 

F = XM(6) + (2.^XM(2)-XM(2)^3«3./4+5./96 .^XM(2)5(X5. )3^SIN(XM(6)) 
1+(1.25D05CXM(2)^^2.-11 ./24 . XXM( 2 ) 3^X4 . +17 ./192 . XXM( 2) xx6 . )XSIN(2.x 
2XM( 

36)) + (13./12 .xXM(2)X5«2.-43./192.XXM(2)X^5. )XSIN(3.XXM(6) ) 

P = XM(1)X(1.-XMC2)XX2.) 

XSP0=-1.5D0xJ(2)x(R/P)^x2.xC0S(XM(3))x(F-XM(6 )+XM(2)xSIN(F)-XMC 2) 
l/2.xsIN(2.xXM(41+F)-SIN(2.X(XM(4)+F))/2.-XM(2)/6 .XSIN(2.XXM(4)+3.X 
2F)) 

RETURN 

END 



XX FUNCTION XSPM x^ 
X5< XX 
XXXXXXXXXXXXXi(XXXXXXXXXXXXXXXX 

xxxxxxxxxxxxxxxxxxxxxxxxxxxxxx 



FUNCTION XSPM(XM) 

DIMENSION XM(6T 
REAL JC4) 

DATA J,R/1 . ,0.00 108276D0, -0.00 000255D0,-0. 0000 0156DO,6 378 .165D0/ 

F = XM(6)+(2.xXMC2)-XM(2)XX3./4+5./96 .XXMC2 )xx5.)xSIN(XM(6)) 

1+(1 .25D0xXM(2)xx2.-ll./24.xXM(2)xx4.+17 ./192 . XXM( 2)xx6 . )X$IN(2.X 
2XMC 

36 ))+( 13./ 12. XXMC 2) XX2.-43./192.XXMC 2) XX5. )XSIN( 3. XXMC6)) 

P = XMCDxci .-XM(2)xx2. ) 

XSPM = -1 .5D0XJ(2)x(R/P)xx2.x$QRT(1.-XM(2)xx 2. )/XM(2)x( (1 .-1.5D0X 
ISINC 

2XM(3))xx2. )X((1.-.25D0XXM(2)XXM(2) )XSIN( F)+XM( 2)/2 . XSINC 2 . XF)+ 
3XMC2) 

4XXM(2)/12 .XSIN(3.XF))+.5D0XSIN(XM(3) )XSIN(XM( 3) )X( (« . 5D0 )x( 1 . + 
51.25D0 

6XXM(2)XXM(2) )xSIN(2.xXM(4)+F)-XM(2)xXM( 2)/8 .xSIN(2.xXM(4)-F)+7 ./6 . 
7x(l ,-XM(2 )xx2./28. )xsIN(2.xxM(4)+3.xF)+.75D0xxM(2)xSIN(2.xxM(4)+4. 
8 X F 

9)+XM(2)xx2 ./8 .xsiN(2.xXM(4)+5.xF)))+9./16 .XJ(2)X(R/P)XX2.X 
1SQRT(1.-XM(2 )xx2.)xSIN(XM(3))XSIN(XM(3))xSIN(2.xXM(4)) 

RETURN 

END 



xxxxxxxxxxxxxxxxxxxxxxxxxxxxxx 
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XX FUNCTION XSPR xx 

XX XX 
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FUNCTION XSPR(XM) 

DIMENSION XM(6) 

REAL J(4) 

DATA J,R/1., 0.00108276 DO, -0. 00 000255D0, -0.0000 0156 DO, 6378. 165D0/ 

F = XM(6)+(2.xXM(2)-XM(2)xx3./4+5./96 .XXM(2)XX5. )XSIN(XM(6)) 
1+(1.25D0XXM(2)XX2.-11./24.XXM(2)XX4.+17 ./192 . XXM( 2) xx6 . )xSIN(2.x 
2XM(6))+(13./12 .xXM(2)xx2.-43./192.XXM(2)XX5.)X3IN(3.xxM(6)) 

P = XM(1)X(1.-XM(2)XX2.) 

RR = P/Cl .+XM(2)XC0S(F)) 

XSPR = -.5D0xJ(2)x(Rxx2./P)X(1 .-1 .5D0XSIN(XM(3))XSIN(XM(3)) )X(1 .+ 
1XMC2) 

2/(l.+SQRT(l .-XM(2)XX2. ))xC0S(F)+2./CSQRT(l .-XM(2)XX2.))XRR/XM(1 
3))+.25XJC2 )xRxx2./PxsINCXM(3))xx2.xC0SC2.xXM(41+2.xF) 
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yi 
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RETURN 

END 



yiyiyiyiyiyiyiyiyiyiyiyiyiyiyiyiyiyiyiyiyiyiyi'Ayiyiyiyiyiyi 

yiyi yiyi 

SUBROUTINE MEAN 

yiyi yiyi 

yiyiyiyiyiyiyiyiyiyiyiyiyiyiyiyiyiyiyiyiyiyiyiyiyiyiyiyiyiyi 
yiyiyiyiyiyiyiyiyiyiyiyiyiyiyiyiyiyiyiyiyiyiyiyiyiyiyiyiyiyi 



SUBROUTINE MEANCXOSC. XM) 



COMPUTES THE MEAN ORBITAL ELEMENTS FROM THE OSCILATORY ELEMENTS 
AND THE FIRST ORDER PERIODIC VARIATIONS BY THE RELATIONS, 

XM = XOSC -XSP(XM) . 

DIMENSION XOSCC6) ,XM(6),XMN(6),XSP(6),TOL(6),DELXM(6) 

DATA TOL/.OIDO, .0000 IDO, .OOOOIDO, . OOOOIDO , . OOOOIDO , .0000 IDO/ 
KOUNT = 0 

yn.., INITIAL GUESS FOR MEAN ELEMENTA. 

DO 10 1=1,6 
10 XM(I) = XOSCCI) 

X... CALCULATE THE SHORT PERIODIC VARIATIONS. 

20 XSPCl) = XSPACXM)- 
XSPC2) = XSPECXM) 

XSPC3) = XSPI(XM) 

XSPC4) = XSPW(XM) 

XSPC5) = XSPO(XM) 

XSPC6) = XSPM(XM) 

KOUNT = KOUNT+1 
DO 50 J=l,6 

XMN(J) = XOSCC J)-XSP(J) 

DELXM(J) = ABS(XMN(J)-XM(J)) 

IF(DELXMCJ) .GT.TOLC J)) XMCJ) = XMN(J) 

50 CONTINUE 
DO 60 K=l,6 

IF(DELXM(K).GT.TOL(K).AND.'KOUNT.LT.30) GO TO 20 
60 CONTINUE 

70 IF(KOUNT.GT.30) WRITE(6,100) DELXM 
90 CONTINUE 

100 F0RMATC/T5, 'CONVERGENCE OF ONE OR MORE MEAN ELEMENTS WAS NOT MET 
l/'THE ABSOLUTE DIFFERENCE IN THE ELEMENTS ARE SHOWN . '//7 E12 . 5) 
RETURN 
END 
yi 
yi 

yi yiyiyiyiyiyi'Ayiyiyiyiyiyiyiyiyiyiyiyiyiyiyiyiyiyiyiyiyiyiyi 

X xxxxxxxxxxxxxxxxxxxxxxxxxxxxxx 

yi XX XX 

^ XX SUBROUTINE GEOP 5CX 

X XX XX 

X XXXXXXXXXXXXXXXXXXXXXXXXXXXXXK 

X xxxxxxxxxxxxxxxxxxxxxxxxxxxxxx 



SUBROUTINE GEOP (XM,DXDT) 

DIMENSION XfK6) , DXDTC6) 

REAL N,NCONS,NCONST,NCON 
REAL JCA) 

DATA J,R/1 ., 0.00 108276 DO, -0. 0000 0255D0, -0.00 000 156 DO, 6 378. 165D0/ 
DATA U/398600./ 

N = SQRT(U/CXM(1)XX3.)) 

P = XM(1)X(1.-XM(2)X^2.) 

DXDT(l) = 0. 

S3 = SIN(XM(3))3(SIN(XM(3)) 

XS = l.-XM(2)5CX2. 

CS = C0S(XM(3)) 

DED = (-3. )/32.XNXJ(2)xJ(2)5C(R/P)xx4.xS3X(14.-15.xS3)5(XM(2)x 

1XS«IN(2.XXMC4)) 
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'A 

A 

A 

A 

A 

A 

A 

A 

A 

A 



DEDO = 3./8 .XN^J(3)5C(R/P;XX3.XSIN(XM(3) )X(4.-5.xS3)xXSxCOS(XM(4)) 
DEDOT = 15./32.XNXJ(4)X(R/P)XX4.XS3X(6 .-7 .XS3)5^(XM(2)XXS) 
1^SIN(2.XXM(4)) 

DXDTC2) = DED-DEDO-DEDOT 

DID = 3./64.5^N).'J(2)xj(2)X(R/P)xx4.^SIN(XM(3)x2. )X(14.-15.XS3)X 
1XM(2)X:«2.KSIN(2.XXM(4)) 

DIDO = 3./S . 5^N3<J(3)X(R/P)XX3.5<C0S(XM(3))X(4.-5.XS3)XXM(2)5<C05(XM(4 
1)) 

DIDOT = 15./64.XNXJ(4)3^(R/P)XX4.XSIN(2.XXM(3))X(6 .-7 .XS3)XXM(2)X 
1XM(2)XSIN(2.5(XM(4)) 

DXDTC3) = DID+DIDO+DIDOT 

DWD = .75D0XN5^J(2)5((R/P)xx2.x(4.-5.5tS3) 

DH = 2.X(14.-15.5(S3)xS3-C2S.-15S.xS3+135.xS3XX2. )XXM(2)XX2. 

DWDO = 3./16 .XN5^J(2)5^J<2.X(R/P)X5(4.X(48-103.XS3 + 215 ./4 . 5^33X5^2 . + 
1(7 .-4.5D05(S3-45./S.)^S3X5(2. )5^XM(2)5(X2.+6 .5^(1 .-1 .5D05fS3)x(4.-5.xS3)x 
2 SQRT(XS)-.25DOXDWXCOS(2.XXM(4))) 

PRINT X, XM(2), XMC3), XM( 4 ) , R , P, N, J ( 3 ) , S3 

DWDOT = 3 ./8 .^NXJ(3)X(R/P)3(X35eC(4 .-5.XS3)x(S3-XMC2)XX2?eCOS(XM(3) 

1) 5(5^2)/(XM(2)XSIN(XM(3)) ) + 2.XSIN(XMC3))X(13.-15.5^S3)XXM(2))5(SIN(XM 
2(4)) 

DHW = S35((6. -7. xS3)-(6. -35. XS3 + 31. 5^33X2^2. )^XM(2)^X2. 

DMDOTT = 15./32.2^N2^J(4)2((R/-P)XK4.2((16.-6 2.3<S5+49.XS3XX2. + .7 5DOX(24.- 
184 .XS3+6 3 . 2 ^S 32 ^^ 2 . )5(Xri(2)3^3<2.+DWWXCOS(2.2^XM(4))) 

DXDT(4) = DWD+DWDO+DWDCT-DNDCTT 
DOD =(-l .5DO)2^N2^J(2)X(R/P)2^2^2 .xCS 

DODO = 1.5D0XN2^J(2)2^?f2.2^(R/P)X2^4.xCS2^(2.25D0 + 1.5xSQRT(XS)-S3 
12((2.5D0+2.25D0 xSQRT( XS ) ) + XM( 2) X5C2 ./ 4 . ( 1 . +1 . 25D0X53 ) -XM( 2 ) 3^X2 . 

2/8 .5€(7 .-15.XS3)3^C0S(2.3^XM(4) )) 

DODOT = 3./8 .XN3^J(3)3^(R/P)X)^3.3C(15.XS3-4)XXM(2)XCS/SIN(XM(3))3^SIN( 
1XM(4)) 

DODOTT = 15./16 .3CN3^J(4)X(R/P)X)^4.3^CS3(( (4.-7.«3)X(l .+1 .5D03^XM(2)3^X2. 

2) -(3.-7 .3^S3)XXM(2)3^3^2.3^COS(XM(4)X2. )) 

DXDT(5) = DOD-DODO-DODOT+DODOTT 

NCON = 3./8.3(N3^J(2)3^)^2. 

NCOMS = N3(J(4)X(R/P)5()^4. 

NCONST = 3./8.3CN3(J(3)3((R/P)3ex3. 

DM = N3C(1 .+1 .5D03^J(2)3^(R/P)XX2.X(1 .-1 .5D0XS3)3^SQRT(XS)) 
DMD=1.5D03^NXJ(2)3^3^2?^(R/P)3^3^43^( (1 . -1 . 5D0XS3 ) 3(X23(XS+( 1 . 25D03(( 1 . - 
12.5D03^ 

133 + 13. /8.XS33^S3) + 5./8. 3^(1. -33+5. /8. 3^333^33 )3(XM( 2 )X3(2+S3/ 16 .3^(14. - 
2 15.3(33)3^(1 .-2. 5D03^Xri( 2)3^X2)3(003 (2. 3(XM(4) ) )3^3QRT(XS) ) 

DMl = (3.-7 .5D03(33+47 ./8 . X33X33 + ( 1 . 5D0-5 , X33+1 17 . /6 . X33X33 ) XXM( 2 ) 3(X2 
1-(1 .+5.X33-101 ./8 .X33X33)/8 .XXM(2)3(X4 . ) 

DM2 = XM( 2) 3(X2. /8.X33X( 7 0. -123. XS3+( 56 . -66 . X33)3(XM( 2)XX2 . )3(C0S( 2 . X 
lXM(4))+27 ./128 .xxM(2)xx4.x33X33xC03(4.XXM(4)) 

DMDO = NC0NX(R/P)xx4./(SQRT(XS))X(3.XDM1+DM2) 

DMDOT = NC0NSTXSIN(XM(3))X(4.-5.XS3)X(1 .-4 .XXM(2)XX2. )/XM(2)X 
15QRT(XS)x 

l^IN(XM(4))-45./128.XNC0N3X(a .-40 .X33+35. X33X33 )xXM(2)xx2.X3QRT(XS) 
DMD0TT=15./64.xnC0NSX33x(6 .-7 .x33)x(2.-5.xXM(2)xx2. )X3QRT(XS)XC0S( 
12.XXM(4)) 

DXDT(6) = DM+DMD+DMDC-DMDOT+DMDOTT 

RETURN 

END 



AAAAAAAAAAAAAAAAAAAAAAAAAAAAAA 
AAAAAAAAAAAAAAAAAAAAAAAAAAAAAA 
AA XX 
XX SUBROUTINE SOLOC XX 
XX XX 
AAAAAAAAAAAAAAAAAAAAAAAAAAAAAA 
AAAAAAAAAAAAAAAAAAAAAAAAAAAAAA 



SUBROUTINE SOLOC ( TVE , T , BLA , DS , AS , AB ) 

DATA EPSl/ 23.5D0 / , PI / 3.14159D0/, TOL/ O.OOOIDO / 
PIl = 0.5D0XPI 
PI2 = 2. X PI 
PIS = 1.5D0X PI 
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RPD = PI / 180 . ODO 

DELT = T - TVE 

EPS = EPSl X RPD 

CLON = 0.985647D0 X RPD X DELT 

IF (CLON.LT. 0.0 ) CLON = PI2 + CLON 

IFCCL0N.GT.PI2) CLON = AMODCCLON, PI2) 

SDS SIN(EPS) SIN ( CLON ) 

DS = ASIN ( SDS ) 

BLAR = BLA ^ RPD 
DO 10 1=1,5 

J = I - 1 

ARG = 0.5D0 J X PI 
TST = ARG - CLON 
TSTl = ABSC TST ) 

IF ( TSTl . LE . TOL ) GO TO 20 
10 CONTINUE 

ETAl = TAN(DS) / TAN(EPS) 

ETA = ASIN(ETAl) 

ETA = ABSC ETA) 

IF (( CLON . LT . PIl ) . AND . ( CLON . GT . 0.0)) AS = ETA 
IFCCCLON.LT.PD.AND.CCLON.GT.PID) AS = PI - ETA 
IF ((CLON. LT. PIS) .AND. (CLON. GT. PI)) AS = PI + ETA 
IF ((CLON. LT.PI2) .AND. (CLON. GT. PIS)) AS = PI2 - ETA 
GO TO 21 

20 AS = CLON 

21 A5 = AS + BLAR 
RETURN 

END 



SUBROUTINE DATPRP xx 

xxxxxxxxxxxxxxxxxxxxxxxxxxxxxx 

xxxxxxxxxxxxxxxxxxxxxxxxxxxxxx 



SUBROUTINE DATPRP ( X,Y,XX,YY ) 

COMMON / PREP / ASUN, DSUN, RAS, DECS 

X = ASUN 

Y = DSUN 

XX = RAS 

YY = DECS 

RETURN 

END 



xxxxxxxxxxxxxxxxxxxxxxxxxxxxxx 

xxxxxxxxxxxxxxxxxxxxxxxxxxxxxx 

XX XX 

XX SUBROUTINE PERIOD xx 
XX XX 

xxxxxxxxxxxxxxxxxxxxxxxxxxxxxx 

xxxxxxxxxxxxxxxxxxxxxxxxxxxxxx 



SUBROUTINE PERIOD ( AO, £S, XI, PA ) 

REAL J2 

COMMON / XXXX / lET 

DATA J2 / 0.00108276D0/, R/ 6S78.165D0/ , U/S986 00 ./ , PI/S . 1 A159D0/ 
PA = 2.)^PI5^(SQRT(A05^A0XA0/U))X(1 ,*-(3.3(J2XRXRX(3.xSIN(XI)xSIN(XI) 

1 -2.0)/(A.xA0XA03<((l .-ES5^ES)XX1 .5D0)))) 

IF ( lET . EQ. 0 ) RETURN 
WRITE (6,100) PA 

100 FCRMAT(/10X, 'INITIAL ANOMALISTIC PERIOD= SF12.S,' SEC',//) 

RETURN 

END 



xxxxxxxxxxxxxxxxxxxxxxxxxxxxxx 
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X 
X 
X 
X 
X 
X 
X 

SUBROUTINE SALT ( D, XM, RM, Z, RE ) 

DIMENSION XMC6) 

DATA ROE/6373. 165D0/, E0E2/6 .693A2E-03/ 

RE = ROE/SQRT(1. + (EOE2/(1.-EOE2))xSIN(D)5(SIN(D)) 

Z = RM + XSPR(XM) - RE 

RETURN 

END 

X 

X 

X xxxxxxxxxxxxxxxxxxxxxxxxxxxxxx 

X xxxxxxxxxxxxxxxxxxxxxxxxxxxxxx 

X XX XX 

XX SUBROUTINE SATLOC xx 

X XX XX 

X xxxxxxxxxxxxxxxxxxxxxxxxxxxxxx 

X xxxxxxxxxxxxxxxxxxxxxxxxxxxxxx 

X 

X 

SUBROUTINE SATLOC (MTYPE, XM, DELT, DSUN, ASUN, DECS , RAS, GMGLT, R) 
DIMENSION XM(6) 

DATA BO/3986 03. 25AD0/,W/. 72921 15E-0 A/. PI/3. 1A159D0/ 

AM = 0.0 
PI02 = 0.5DOXPI 
PI2 = 2.XPI 
P2 = PI2 
PI03 = 3.XPI/2. 

AG = AMOD ( W5^DELT , PI2 ) 

A = XM(l) 

CE = XMC2) 

Cl = XM(3) 

G = AMOD ( XM(A) , PI2 ) 

H = AMOD ( XMC51 , PI2 ) 

CALL POSVEL(A,CE,CI,AM,G,H,BO,R1,R2,R3,V1,V2,V3,R,V) 

DECS = ATAN ( R3/SQRTC RlxRl + R2^R2 )) 

RAS = ATAN2 ( R2 , R1 ) 

IF -( MTYPE .EQ. 0 ) GO TO 10 

ALP = TAN(DSUN) x COS(CI) / SIN(CI) 

ALP = ASIN(ALP) 

ALP = ABS(ALP) 

DECS = DSUN 

IF ( MTYPE . EQ . 1 ) DECS = - DSUN 
X U1 = ASINC SIN(DSUN) / SINCXMC3))) 

X U2 = ASIN ( SIN(DECS)/SIN(XMC3))) 

T1 = ASUN + PI02 
T2 = ASUN - PI02 
Wil = G 
OMA = H 

IF (T2 .LT. 0.0 ) GO TO 20 
- IF ( T1 .GT. PI2 ) GO TO 30 

IF (COMA .LE. Tl) .AND. (OMA .GE. T2)) GO TO AO 
GO TO 50 

20 T3 = PI2 + T2 

IF (((OriA.LE.Tn .AND. (OMA.GE.O. 0)) . OR . ( ( OMA . L E . P2) . AND . ( OMA . GE . T3 
1 ))) GO TO AO 

GO TO 50 

30 T3 = Tl - PI2 

IF(( (0MA.GE.T2) . AND .( OMA . LE . PI2) ). OR .(( OMA . GE . 0 .). AND .( OMA . LE . T3 
1 ))) GO TO AO 

GO TO 50 

AO IF ( Cl . GE . PI02 ) GO TO A1 
RAS = OMA + ALP 

IF ( DECS . LT . 0.0 ) RAS = OMA - ALP 
GO TO A2 

A1 RAS = OMA - ALP 



xxxxxxxxxxxxxxxxxxxxxxxxxxxxxx 

XX XX 

XX SUBROUTINE SALT xx 

XX 5^x 

xxxxxxxxxxxxxxxxxxxxxxxxxxxxxx 

xxxxxxxxxxxxxxxxxxxxxxxxxxxxxx 
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'A 

X 

X 

X 

n 

yi 



yi 

yi 

yi 

yi 



y^ 

X 

yi 

yi 



yi 

yi 





IF ( 


DECS 


. LT . 


0.0 


) RAS 


= OMA 


+ ALP 


42 


IF ( 


MTYPE 


. EQ 


. 1 : 


1 RAS = 


RAS + PI 




RAS = 


AMOD 


( RAS 


, PI2 ) 








GO TO 


10 












50 


IF ( 


Cl . 


GE . PI02 ) GO TO 


51 






RAS = 


OMA 


+ PI - 


ALP 










IF (DECS . 


LT . 


0.0 : 


1 RAS = 


OMA + PI + 




GO TO 


42 












51 


RAS = 


OMA 


+ PI + 


ALP 










IF ( 


DECS 


. LT . 


0.0 


) RAS 


= OMA 


+ PI 




GO TO 


42 












10 


ELNG 


= RAS 


- AG 












IF (ELNG . 


LT. 0. 


0 ) ELNG = 


PI2 + 


ELNG 



ELNG = AM0D(ELNG,PI2) 

GMGLT = ASIN ( 0 . 9792D0?(SIN( DECS) + 0 . 2028D05(C0S( DECS ))^COS( ELNG - 
1 (291.5(PI/180.))) 

RETURN 

END 



yiyiyiyiyiyiyiyiyiyiyiyiyiyiyi^yi^iyiyiyiyi^yi^yiyiyiyiyi 
yiyiyiyiyiyiyiyiyiyiyi^yiyiyiyiyiyiyiyiyiyiyiyiyiyiyiyiyiyi 
yiyi yiyi 

XX SUBROUTINE DRAG xx 

XX XX 

xxxxxxxxxxxxxxxxxxxxxxxxxxxxxx 

xxxxxxxxxxxxxxxxxxxxxxxxxxxxxx 



SUBROUTINE DRAGC XM, DEL, TT, DELMG ) 

REAL N, 10, II, 12, 13, 14, 15, 16, 17 

COMMON / SOLAR / F107 , FBAR, AKP, TVE 

COMMON / PREP / ASUN,DSUN,RAS,DECS 

DIMENSION XM(6), SUNC4) , SATC2) , GE0C3) , BF(8) 

DATA BO/ 3986 03. 254D0/,BLA/30./, PI/3. 14159 DO/, EOE/ .00335D0/ 

DATA IYR/1981D0/,TVD/1 ./ , THR/ 17 ,/ , TM/14 ./ , TS/56 . 243D0/ 

PI2 = 2. X PI 

SUNC3) = 23.5D0XCPI/180.) 

SUNC4) = 1.0 

GEO(l) = F107 

GE0C2) = FBAR 

GE0C3) = AKP 

N = SQRT(B0/(XM(1)XX3)) 

.TCOR = AMOD ( XM(6) , PI2 ) 

TCOR = TCOR / N 

TTT = TT - ( TCOR / 86400. ) 

TTR = 1721044. + 367XIYR - (7XIYR/4) + TVD + (( THRX3600. + TMX60. 
1 + TS ) / 86400. ) - 2400001. ODO 

DELT = TTT - TTR 

CALL SOLOC(TVE,TTT,BLA,DSUN,ASUN,ABLG) 

MTYPE = 0 

CALL SATLOC(MTYPE,XM, DELT, DSUN,ASUN, DECS ,‘RAS, GMGLT, RM) 

STR = XM(6) 

XM(6) = 0.00 

CALL SALTC DECS , XM, RM, HP, RE ) 

XM(6) = STR 

CALL DATPRPC SUN( 1) , SUN( 2) , SATC 1) , SAT( 2) ) 

CALL ATMDENC TTT, SUrJ, SAT, GMGLT, HP, GEO, RHO,DRDH,SH) 

RHOP = RHO X ( 10. XX9. ) 

SH? = SH X O.OOIDO 

CALL JAC60(DECS,DSUN,RAS,ASUN,BLA,HP, F107,RH0P,SHP) 

MTYPE = 1 

CALL SATLOC(MTYPE,XM, DELT, DSUN,ASUN, DECS, RAS, GMGLT, RM) 

CALL DATPRP ( SUN( 1 ) , SUN ( 2 ) , S AT( 1 ) , SATC 2 ) ) 

CALL ATMDENC TTT, SUN, SAT, GMGLT, HP, GEO, RHO, DRDH,SH) 

RHOMIN = RHO X ( 10.XX9 ) 

SHMAX = SH X O.OOIDO 

CALL JAC6 0(DECS,DSUN,RAS,ASUN, BLA, HP, FI 07, RHOMIN, SHMAX) 

MTYPE = 2 

CALL SATLOC(MTYPE,XM,DELT,DSUN,ASUN,DECS,RAS,GMGLT,RM) 

CALL DATPRP ( SUNCl), SUNC2), SAT(l), SATC2) ) 

CALL ATMDEN(TTT,SUN,SAT,GMGLT,HP,GEO,RHO,DRDH,SH ) 
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. . . 



X 

X 

X 

X 



RHOMAX=RHO yi C10.X3(9 ) 

SHMIN = SH^^O.OOIDO 

CALL JAC60CDECS, DSUN, RAS, ASUN, BLA, HP, F107 , RHOMAX, SHMIN) 

A = SIN(DSUN)3€SIN(XM(3))9<SIN(XM(4))+COS(DECS)3((COS(XM(5)-ABLG)X 
1 COS(XM(4))-COS(XM(3))9(SIH(XM(5)-ABLG)3(SIN(XM(4) )) 

B = SIN(DSUN)XSIN(XM(3) )XCOS(XM(4))-COSCDECS)3^(COS(XM(5)-ABLG)B( 

1 SIN(XM(4))+C0S(XM(3) )^SIN(XM(5)-ABLG)5^COS(XM(4) )) 

F = (RHOMAX - RHOMIN) / (RHOMAX + RHOMIN ) 

RHOO = RHOP - 0 .5D03^(RHOMAX + RHOMIN) X F X A 

Z = XM(1) yi XMC2) / SHP 

Z1 = -!.)-( F3(A / (1. + F X A )) 

Z2 = ( 1. / ( 1. + F 5^ A )) - Z 5^ ( 1. - 0.5D0^Z) - 0.5D0 
TMl = -Zly^Zl - 2. Z^Z^ZZ 
CTHAVG =0.0 

IF ( TMl .GT. 0.00 ) CTHAVG=(Z1+SQRT(TM1))/(Z3^Z) 

SHNM = ABS( SHMAX - SHMIN ) 

H-AVG = 0.5D03^(SHMAX+SHMIN)+ ( SHNM / ( SHMAX+SHMIN) )XA^CTHAVG 
BETA = 1.0 / HAVG 

C = 0.5DOXBETA3^EOE3^SIN(XMC3))3^SIN(XM(3))X(RE+HP) 

CC = C c 

ZZ - BETA X XM(1) X XM(2) 

ALP = 0.00 
MO = 1 
NN = 3 

CALL MMBSIR(ZZ,ALP,NN,MO,BF,NZ) 

10 = BF(1) 

11 = BF(2) 

12 = BF(3) 

13 = BF(4) 

14 = BF(5) 

15 = BFC6) 

16 = BFC7) 

17 = BF(8) 

FA = FXA 
FB = FXB 

E = XM(2) 

C2W = C0SC2.XXMC4)) 

S2W = SIN(2.XXM(4)) 

C4W = COS(4.XXM(4)) 

S4W = SIN(4.XXM(4)) 

DELMG=-DELX(SQRT(BOxXM(1)))xRHOOxEXP(-ZZ-CxC2W)x(( 1 .+0.25D0X 
lCxC)x 

1(I0+EXI1)+FAX((1 .+0.25D0XCC)x(I1+ExI2) )+0 . 5D0XCX( (2 . XI2-EX( I 1-3 . x 
2I3))+FAX( I l+I 3+2. XEXI4))XC2W-0.5D0XCXFDX(( 11-13 )+2,xEX( 12-14) )x 

3 S2W +0 .125D0XCXCX((2.XI4-£X(3.XI5-5.XI5))+FAX((I3+I5)+EX(3.XI6-I2 

4 )))XC4W + 0.125D0XCXCXFBX((I5-I3)+EX(I2-4.XI4+3.XI6))XS4W) 

XTl = DELX(SQRT(BOXXM(1 )))xRHOO 

XT2 = DELMG / XTl 

XT3 = EXP(-ZZ)x(I0 + EXIl) 

XT4 = XT2 / XT3 

DELMG = 1000. X DELMG / 9.8D0 

PRINT X, ZZ, RHOO, XTl, XT2, XT3, XT4, DELMG , DEL 

RETURN 

END 



XX5(^9(XXXXXXXXXXXXXXXXXXXXXXXXX 
B(XXXXXXXXXXXXX3€XXXXX3€XX^XXXXXX 
XX XX 

XX SUBROUTINE JAC60 xx 

XX XX 

xxxxxxxxxxxxxxxxxxxxxxxxxxxxxx 

xxxxxxxxxxxxxxxxxxxxxxxxxxxxxx 



SUBROUTINE JAC60 ( D, DS , A , AS, BLA, Z, FI 07 , RHO , SH ) 

BL = BLAX3. 14159D0 / 180. 

XP = -16.021D0 - 0.001985D0XZ + 6.363D0xEXP( -0.0026D0XZ) 
RO = 10. xxXP 

CSPSI =’SIN(D)XSIN(DS) + COS( D)XC0S( DS)XC0S( A-AS-BL ) 

F = (O.SDOXd. + CSPSI))XX3. 

RHOl = 1. + 0.19D0X(EXP(0.0055D0XZ) - 1.9D0) XF 
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RHO = O.OIDO F107 X RO X RHOl ^ (10.3^^12) 

DRODZ = 2.3026D03CR03((-.001985D0 - . 016 544D03^EXP( 0026 D03(2) ) 
DR1D2 = .0010^*5D03(FXEXP( .0055D0XZ) 

SH = “(ROXRH01V(R03(DR1DZ+RH013(DRODZ) 

RETURN 

END 



SUBROUTINE THRSTC XM, DXM, IRC, K, PA, DWOA ) 

DIMENSION XM(6) , DXMC6), IRV(IO), DA(IO) 

COMMON / OADJ / NOA, IRV, DA, SPI, WT, PSENS 
DATA BO / 398603. 254D0/, PI/3 . 1 A159D0/ 

DO 10 I = 1,6 
DXM(I) = 0.0 
10 CONTINUE 
DWOA =0.0 

IF ( IRC .EQ. 0) RETURN 

THT = XM(6) + 2.3^XM(2)XSIN(XM(6)) + ( 5 ./4 . )XXM( 2) 3(XM( 2) 3^SIN( 2 . XXM 
1 ( 6 )) 

R = XM(1)X(1. - XM(2)XXM(2))/(1 . + XM(2)3(COS(THT) ) 

V2 = BO ((2./R) - (l./XMCl))) 

V = SQRTCV2) 

DVI = DACK)3(B0/(2.XXM(1)3(XM(1)3(V) 

DXM(1)= DA(K) 

DXMC2) = 2.3(((COSCTHT)+XM(2))/V)3(DVI 
DXM(A) = 2.XSIN(THT)XDVI/(XM(2)3CV) 

DXM(6) = (2.3^PI/PA)-(1./(XM(1)X(SQRT(1 .-XM(2:3^XM(2)))3(V))3(2.3( 

1 SIN(THT)3C((XM(1)X(1.-XM(2)3(XM(2))/XM(2)) + RXXM(2))XDVI 

PA = 3.3CPAXXM(1)3C(V/B0)3^DVI 
DWOA = (WT/(SPl3(9.8D0))3a000 .XDVI 
. IF (PSEMS .EQ. 0) THEN 

PRINT THT,R,V,IRV(K),DA(K),DVI,DXM(2),DXM(A),DXM(6),PA, 



X 

X 

X 

X 

X 

X 

X 

X 

X 

X 

X 



xxxxxxxxxxxxxxxxxxxxxxxxxxxxxx 

xxxxxxxxxxxxxxxxxxxxxxxxxxxxxx 



XX 

XX 

XX 



xxxxxxxxxxxxxxxxxxxxxxxxxxxxxx 

xxxxxxxxxxxxxxxxxxxxxxxxxxxxxx 



SUBROUTINE THRST 



XX 

XX 

XX 



DWOA 



END IF 
RETURN 
E!1D 
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APPENDIX B 



SAMPLE INPUT AND OUTPUT 



05000,001,0,1,0,025 

2.2, .0000500, 1.0, 230. ,100. ,09100. 

ICO., 100., 2.0 ,43222.7382 

1 10947 U 90070 78 060 A 0 

2 10947 44619.98716775 150.0000 

3 13947 15.71356734 .000783912 

4 10947 000000000-0 0 

5 10947 01.06000000 -701.8652 E-00 



1.0 014791 1015 0 0 ++ 

095.0000 200.0000 .0100000 095.0000 ++ 



1 , 



00112 

0 



-3.81308 -0.23959E-00 
0 0 



++ 

++ 

++ 



cccc 


X.X 


CD 


1.0 


A 


0.0 


D 


0.0 


SPI 


0.0 


N 


0.0 


WT 


0.0 


F107 


0.0 


FBAR 


0.0 


AKP 


0.0 


TVE 


0.0 


UJD 


0.0 


ETU 


0.0 


HO 


0.0 


GO 


0.0 


ES 


0.0 


XI 


0.0 


RND 


0.0 


ESD 


0.0 


AO 


0.0 


ADOT 


0.0 



XXXXX) 



axxxxxx 



000001.000000000 

000000.000010000 

000000.100000000 

000200.000000000 

000005.000000000 

000100.000000000 
000050.000000000 

000050.000000000 
000000.000000000 

038422.000000000 

043223.000000000 
000000 . 000000000 
000000. ooooocooo 
000000 . 000000000 
000000.010000000 
000022.500000000 
000000.000010000 
000000.000000000 
000001 . 025940000 
000000.000000000 



1000 , 0.1 
2000, 0.3 
3000, 0.5 



xxxxxx.xxxxxxxxx 

000003.500000000 
000000 . 001000000 
000010.000000000 

003000.000000000 

000200 . 000000000 
100000 . 000000000 
000300.000000000 

000300.000000000 

000007.000000000 

039422.000000000 

044223.000000000 

000179.000000000 

000360.000000000 
000360 . 000000000 
000000 .100000000 

000157.500000000 
000000.000000000 
000000 . 000000000 
000001 . 168470000 
000000 . 000000000 



yi 

^ ONLY THE ABOVE IS ACTUAL INPUT FILE DATA. AN INPUT ELEMENT 

5^ DESCRIPTION IS PRESENTED NEXT TO AID THE USER IN PLACING VALUES x 

yi 

INPUT FILE FOR PROPELLANT LONGEVITY MODEL (VARIABLE ORGANIZATION) 

lETRl, IETR2, NOA, PSENS, lET, PSIZE 
CD, A, D, SPI, W, WT, 

F107, FBAR, AKP, TVE, 

1 SKIP LINE ONE NAVSPASUR FIVE CARD DATA 

2 8X, UJD, IX, ETU, IX, HO, IX, GO, IX, ES, IX, XI 

3 20X, RND, 19X, ESD 

4 SKIP LINE FOUR NAVSPASUR FIVE CARD DATA 

5 3X, AC, IX, ADOT 

1 

NOTE: WA.RNING OMIT ALL ORBIT ADJUST DATA IF 
NO ORBIT ADJUSTS ARE SCHEDULED CNOA=0) 



IRVCl), 


DA(l) 


IRVC2), 


DAC2) 


IRV(3), 


DA(3) 


IRVCK), 


DAU) 



I 



SENSITIVITY ANALYSIS TABLE 



CHAR SNSAR 
A4 F3.1 



SNSMIN 
F16 .8 



SNSMAX 

F16.8 
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INITIAL INPUT CONDITIONS 



ATMOSPHERIC INPUT PARAMETERS 

D (RELATIVE ATMOSPHERIC ROTATION)* = 1.0000 

F107 (DECIMETRIC SOLAR FLUX) = 100.00 

FBAR (90 DAY AVG . F107) = 100.00 

AKP (GEOMAGNETIC INDEX) = 2.000 

TVE (TIME OF VERNAL EQUINOX PASSAGE) = A3222 .733200000 MODIFIED JULIAN DAYS 



SATELLITE PHYSICAL PARAMETERS 

CD (DRAG COEFFICIENT) = 2.20 UNITLESS 

A (CROSSECTIONAL AREA) = 0.000050000 SQUARE KILOMETERS 
SPI (MOTOR SPECIFIC IMPULSE) = 230.00 SECONDS 

W (INITIAL FUEL MASS) = 10.00 KILOGRAMS 

WT (INITIAL SATELLITE MASS) = 9100.00 KILOGRAMS 



NAVSPASUR ORBITAL ELEMENTS 

UJD (EPOCH) = AA619. 98716775 MEAN JULIAN DAYS 
ETU (MEAN ANOMALY) = 150.0000 DEGREES 

HO (MEAN RIGHT ASCENSION OF THE ASCENDING NODE) = 95.0000 DEGREES 

GO (MEAN ARGUMENT OF PERIGEE) = 200.0000 DEGREES 
XI (MEAN INCLINATION) = 95.0000 DEGREES 

RND (RATE OF CHANGE OF MEAN MOTION) = 0.000783912 REVOLUTTIONS PER DAY SQUARED 

ES (MEAN ECCENTRICITY) = 0.0100 UNITLESS 

ESD (ECCENTRICITY DECAY RATE) = -0.239590000 1/SECONDS 

AO (KAULA SEMI-MAJOR AXIS) = 1,06000 MEAN EARTH RADII 

ADOT (SEMI-MAJOR AXIS DECAY RATE) = -701.865200000 



CONTROL INPUT PARAMETERS 

MAXIMUM REVOLUTIONS PER ITERATION 5000 

NUMBER OF SCHEDULED ORBIT ADJUSTS TO SEMI-MAJOR AXIS 0 



SENSITIVITY ANALYSIS FOR 
LOW EARTH ORBIT SATELLITES 
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TABLE OF INPUT PARAMETERS SCANNED FROM MINIMUM TO 
MAXIMUM WITH THE INDICATED INCREMENT PER RUN. 

THE GRANULARITY IS 1/ 25. 

ALL OTHER INPUT PARAMETERS ARE RETURNED TO THEIR 
INITIAL VALUES AT EACH INCREMENT. 



INPUT PARAMETER MINIMUM 



MAXIMUM 



INCREMENT 



CD 



1.000000000 



3.500000000 



0.100000000 



SENSITIVITY SCAN RESULTS 



REV NUMBER 


TIME(DAYS) 


495 


31.58 


448 


28.57 


409 


26.08 


376 


23.97 


348 


22.18 


324 


20.65 


303 


19.30 


285 


18.15 


268 


17.07 


254 


16.17 


241 


15.54 


229 


14.57 


218 


13.87 


208 


13.23 


200 


12.72 


191 


12.14 


184 


11.70 


177 


11.25 


170 


10.30 


154 


10.42 


159 


10.10 


154 


9.78 


149 


9.46 


144 


9.14 


140 


8.88 


136 


8.63 



FOR ABOVE SPECIFICS 

FRACTION OF SCAN 

1 .000000000 

1.100000000 

1.200000000 

1.300000000 
>.^00000000 

1.500000000 
1.600000000 

1.700000000 
1.800000000 
1.900000000 

2.000000000 

2.100000000 

2.200000000 

2.300000000 

2.400000000 

2.500000000 
2.600000000 

2.700000000 
2.800000000 
2.900000000 

3.000000000 

3.100000000 

3.200000000 

3.300000000 

3.400000000 

3.500000000 
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05000.001.0. 1.0.025 

2.2, .0000500, 1.0, 250. ,001. ,09100 

100., 100., 2.0 ,43222.7382 

1 10947 U 90070 78 060 A 0 

2 10947 44619.98716775 150.0000 095. 

3 10947 15.71356754 .000733912 1.0 

4 10947 000000000-0 0 

5 10947 01.06000000 -701.8652 E-00 



cccc 


X.X 


CD 


0.0 


A 


0.0 


D 


0.0 


SPI 


0.0 


W 


0.0 


WT 


0.0 


F107 


0.0 


FBAR 


0.0 


AKP 


0.0 


TVE 


0.0 


UJD 


0.0 


ETU 


0.0 


HO 


0.0 


GO 


0.0 


ES 


0.0 


XI 


0.0 


RND 


0.0 


ESD 


0.0 


AO 


1.0 


ADOT 


0.0 



xxxxxx . xxxxxxxxx xxxxxx . 

000001.000000000 000003. 

000000.000010000 000000 . 

000000.100000000 000010 . 

000200.000000000 003000. 

000005.000000000 000200. 

000100.000000000 100000 . 

000050.000000000 000300. 

000050.000000000 000300 . 

000000.000000000 00D007 . 

038422.000000000 039422. 

043223.000000000 044223. 

000000.000000000 000179 . 

000000.000000000 000360. 

000000.000000000 , 000360. 

000000.010000000 000000 . 

000022.500000000 000157 . 

000000.000010000 000000 . 

000000.000000000 000000 . 

000001.025940000 000001. 

000000.000000000 000000 . 



1.0 014791 1015 0 0 ++ 

0000 200.0000 .0100000 095.0000 ++ 
0112 -3.81308 -0.25959E-00 0 ++ 

0 0 0 ++ 

++ 

xxxxxxxxx 

500000000 

001000000 

000000000 

000000000 

000000000 

000000000 

000000000 

000000000 

000000000 

000000000 

000000000 

000000000 

000000000 

000000000 

100000000 

500000000 

000000000 

000000000 

168470000 

000000000 



1000 , 0.1 
2000, 0.3 
3000, 0.5 



X 

X ONLY THE ABOVE IS ACTUAL INPUT FILE DATA. AN INPUT ELEMENT 
X DESCRIPTION IS PRESENTED NEXT TO AID THE USER IN PLACING VALUES 
X 

xxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxx 



INPUT FILE FOR PROPELLANT LONGEVITY MODEL (VARIABLE ORGANIZATION) 

lETRl, IETR2, NOA, PSENS, lET, PSI2E 
CD, A, D, SPI, W, WT, 

F107, FBAR, AKP, TVE, 

1 SKIP LINE ONE NAVSPASUR FIVE CARD DATA 

2 3X, UJD, IX, ETU, IX, HO, IX, GO, IX, ES, IX, XI 

3 20X, RND, 19X, ESD 

4 SKIP LINE FOUR NAVSPASUR FIVE CARD DATA 

5 8X, AO, IX, ADOT 
IRVCl), DA(1) I 

IRVC2), DA(2) I NOTE: WARNING OMIT ALL ORBIT ADJUST DATA IF 
IRVC3), DA(3) I NO ORBIT ADJUSTS ARE SCHEDULED (N0A=0) 



IRV(K), DA(K) I 

SENSITIVITY ANALYSIS TABLE 

CHAR SNSAR SNSMIN SNSMAX 

A4 F3.1 F16.8 F16.8 

xxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxx 
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X >K X 



INITIAL INPUT CONDITIONS 



ATMOSPHERIC INPUT PARAMETERS 

D (RELATIVE ATMOSPHERIC ROTATION) = 1.0000 

F107 (DECIMETRIC SOLAR FLUX) = 100.00 

F2AR (90 DAY AVG . F107) = 100.00 

AKP (GEOMAGNETIC INDEX) = 2.000 

TVE (TIME OF VERNAL EQUINOX PASSAGE) = ^3222.738200000 MODIFIED JULIAN DAYS 



SATELLITE PHYSICAL PARAMETERS 

CD (DRAG COEFFICIENT) = 2.20 UNITLESS 

A (CROSSECTIONAL AREA) = 0.000050000 SQUARE KILOMETERS 
SPI (MOTOR SPECIFIC IMPULSE) = 230.00 SECONDS 

W (INITIAL FUEL MASS) = 1.00 KILOGRAMS 

HT (INITIAL SATELLITE MASS) = 9100.00 KILOGRAMS 



NAVSPASUR ORBITAL ELEMENTS 

UJD (EPOCH) = 4A619. 98716775 MEAN JULIAN DAYS 
ETU (MEAN ANOMALY) = 150.0000 DEGREES 

HO (MEAN RIGHT ASCENSION OF THE ASCENDING NODE) = 95.0000 DEGREES 

GO (MEAN ARGUMENT OF PERIGEE) = 200.0000 DEGREES 
XI (MEAN INCLINATION) = 95.0000 DEGREES 

RMD (RATE OF CHANGE OF MEAN MOTION) = 0.000783912 REVOLUTTIONS PER DAY SQUARED 

ES (MEAN ECCENTRICITY) = 0.0100 UNITLESS 

ESD (ECCENTRICITY DECAY RATE) = -0.239590000 l/SECONDS 

AO (KAULA SEMI-MAJOR AXIS) = 1.06000 MEAN EARTH RADII 

ADOT (SEMI-MAJOR AXIS DECAY RATE) = -701.865200000 



CONTROL INPUT PARAMETERS 

MAXIMUM REVOLUTIONS PER ITERATION 5000 

NUMBER OF SCHEDULED ORBIT ADJUSTS TO SEMI-MAJOR AXIS 0 



SENSITIVITY ANALYSIS FOR 
LOM EARTH ORBIT SATELLITES 
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TABLE OF INPUT PARAMETERS SCANNED FROM MINIMUM TO 
MAXIMUM HITH THE INDICATED INCREMENT PER RUN. 

THE GRANULARITY IS 1/ 25. 

ALL OTHER INPUT PARAMETERS ARE RETURNED TO THEIR 
INITIAL VALUES AT EACH INCREMENT. 



INPUT PARAMETER MINIMUM MAXIMUM 



INCREMENT 



AO 



1.0259A0000 



.16SA70000 0.005701200 



SENSITIVITY SCAN RESULTS 



MBER 


TIMECDAYS) 


1 


0.00 


1 


0.00 


1 


0.00 


3 


0,12 


6 


0,31 


12 


0.70 


22 


1.34 


40 


2.51 


69 


4.42 


115 


7.47 


186 


12.21 


298 


19.76 


47 3 


31.65 


747 


' 50.42 


1149 


78.19 


1652 


113.33 


2033 


140.55 


2270 


158.15 


2514 


176.50 


2942 


208.13 


3349 


238.73 


3643 


262.00 


4026 


291.32 


4373 


318.80 


4555 


334.53 


4712 


348.62 



FOR ABOVE SPECIFICS 

FRACTION OF SCAN 
1.0259^0000 
1.0316A1200 
1 .0373A2A00 
1 .0A30A3600 
1 .0A87AAS00 
1 .05AAA6000 
1.0601A7200 
1 .0658A8A00 
1 .0715A9600 
1.077250800 
1.082952000 
1.088653200 
1 .09435AA00 
1.100055600 
1 .105756300 
1.111458000 
1.117159200 
1 . 122860400 
1.1285616C0 
1.134262300 
1.139964000 
1.145665200 
1.151366400 
1.157067600 
1.162768300 
1.163470000 
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05000, 
2 . 2 , . 
100 ., 

1 10947 

2 10947 

3 10947 

4 10947 

5 10947 
CCCC X. 
CD 0. 



001.0. 1.0.025 

0000500, 1.0, 230. ,100. ,09100. 

100., 2.0 ,43222.7382 

U 90070 78 060 A 0 1.0 014791 1015 0 0 ++ 

44619.93716775 150.0000 095.0000 200.0000 .0100000 095.0000 ++ 



15.71356734 .000783912 1.00112 
000000000-0 0 0 
01.06000000 -701.8652 E-00 



-3.81308 -0.23959E-00 
0 0 



++ 

++ 

++ 



X xxxxxx.xxxxxxxxx 
0 000001.000000000 
A 0.0 000000.000010000 
D 0.0 000000.100000000 
SPI 0.0 000200.000000000 
W 0.0 000005.000000000 
WT 0.0 000100.000000000 
F107 1.0 000050.000000000 
FBAR 0.0 000050.000000000 
AKP 0.0 000000. 000000000 
TVE 0.0 038422. OOOOOOOGO 
UJD 0.0 043223.000000000 
ETU 0.0 000000.000000000 
HO 0.0 000000.000000000 
GO 0.0 000000.000000000 
ES 0.0 000000.010000000 
XI 0.0 000022.500000000 
RNO 0.0 000000.000010000 
ESD 0.0 000000.000000000 
AO 0.0 000001.025940000 
ADOT 0.0 000000.000000000 



XXXXXX.XXXXXXXXX 

000003.500000000 
000000.001000000 

000010.000000000 

003000.000000000 

000200.000000000 

100000.000000000 

000300.000000000 

000300.000000000 

000007.000000000 

039.422.000000000 

044223.000000000 

000179.000000000 

000360.000000000 

000360.000000000 

000000.100000000 

000157 .500000000 
000000.000000000 
000000.000000000 
000001.168470000 
000000.000000000 



1000 , 0.1 
2000, 0.3 
3000, 0.5 



X 

X ONLY THE ABOVE IS ACTUAL INPUT FILE DATA. AN INPUT ELEMENT 
X DESCRIPTION IS PRESENTED NEXT TO AID THE USER IN PLACING VALUES 
X 

xxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxx 

INPUT FILE FOR PROPELLANT LONGEVITY MODEL (VARIABLE ORGANIZATION) 

lETRl, IETR2, NOA, PSENS, lET, PSIZE 
CD, A, D, SPI, W, WT, 

F107, FBAR, AKP, TVE, 

1 SKIP LINE ONE NAVSPASUR FIVE CARD DATA 

2 3X, UJD, IX, ETU, IX, HO, IX, GO, IX, ES, IX, XI 

3 20X, RND, 19X, ESD 

4 SKIP LINE FOUR NAVSPASUR FIVE CARD DATA 

5 SX, AO, IX, ADOT 
IRV(l), DACl) I 

IRVC2), DA(2) I NOTE: WARNING OMIT ALL ORBIT ADJUST DATA IF 
IRVC3), DA(3) I NO ORBIT ADJUSTS ARE SCHEDULED (NOA=0) 



IRV(K), DACK) I 

SENSITIVITY ANALYSIS TABLE 

CHAR SNSAR SNSMIN SNSMAX 

A4 F3.1 F16.8 F16.8 

XXXXXXXXXXXXXXXXXXXXX^^XXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXX 
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X X X 



INITIAL INPUT CONDITIONS 



ATMOSPHERIC INPUT PARAMETERS 

D (RELATIVE ATMOSPHERIC ROTATION) = 1.0000 

F107 (DECIMETRIC SOLAR FLUX) = 100.00 

FBAR (90 DAY AVG . F107) = 100.00 

AKP (GEOMAGNETIC INDEX) = 2.000 

TVE (TIME OF VERNAL EQUINOX PASSAGE) = A3222 . 738200000 MODIFIED JULIAN DAYS 



SATELLITE PHYSICAL PARAMETERS 

CD (DRAG COEFFICIENT) = 2.50 UNITLESS 

A (CROSSECTIONAL AREA) = 0.000050000 SQUARE KILOMETERS 
SPI (MOTOR SPECIFIC IMPULSE) = 230.00 SECONDS 

W (INITIAL FUEL MASS) = 10.00 KILOGRAMS 

WT (INITIAL SATELLITE MASS) = 9100.00 KILOGRAMS 



NAVSPASUR ORBITAL ELEMENTS 

UJD (EPOCH) = AA619. 98716775 MEAN JULIAN DAYS 
ETU (MEAN ANOMALY) = 150.0000 DEGREES 

HO (MEAN RIGHT ASCENSION OF THE ASCENDING NODE) = 95.0000 DEGREES 

GO (MEAN ARGUMENT OF PERIGEE) = 200.0000 DEGREES 
XI (MEAN INCLINATION) = 95.0000 DEGREES 

RND (RATE OF CHANGE OF MEAN MOTION) = 0.000783912 REVOLUTTIONS PER DAY SQUARED 

ES (MEAN ECCENTRICITY) = 0.0100 UNITLESS 

ESD (ECCENTRICITY DECAY RATE) = -0.239590000 1/SECONDS 

AO (KAULA SEMI-MAJOR AXIS) = 1.06000 MEAN EARTH RADII 

ADOT (SEMI-MAJOR AXIS DECAY RATE) = -7 01 . 865200000 



CONTROL INPUT PARAMETERS 

MAXIMUM REVOLUTIONS PER ITERATION 5000 

NUMBER OF SCHEDULED ORBIT ADJUSTS TO SEMI-MAJOR AXIS 0 



SENSITIVITY ANALYSIS FOR 
LON EARTH ORBIT SATELLITES 
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TABLE OF INPUT PARAMETERS SCANNED FROM MINIMUM TO 
MAXIMUM WITH THE INDICATED INCREMENT PER RUN. 

THE GRANULARITY IS 1/ 25. 

ALL OTHER INPUT PARAMETERS ARE RETURNED TO THEIR 
INITIAL VALUES AT EACH INCREMENT. 



INPUT PARAMETER 
F107 



MINIMUM 

50.000000000 



MAXIMUM 

300.000000000 



INCREMENT 

10.000000000 



SENSITIVITY SCAN RESULTS 



REV NUMBER 


TIMECDAYS) 


392 


24.99 


324 


20.65 


276 


17 .58 


241 


15.34 


213 


13.55 


191 


12.14 


174 


11.06 


159 


10.10 


146 


9.27 


136 


8.63 


127 


8.05 


119 


7.54 


112 


7.10 


105 


6.65 


100 


6.33 


95 


6.01 


90 


5.69 


86 


5.43 


82 


5.13 


79 


4.99 


76 


4.79 


73 


4.60 


70 


4.41 


67 


4.22 


65 


4.09 


63 


3.96 



FOR ABOVE SPECIFICS 

FRACTION OF SCAN 

50.000000000 

60.000000000 

70.000000000 

80.000000000 

90.000000000 

100.000000000 

110.000000000 

120.000000000 

130.000000000 
1 ^ 0 . 000000000 

150.000000000 

160.000000000 

170.000000000 

180 .000000000 

190.000000000 

200.000000000 

210.000000000 

220 .000000000 

230.000000000 
2 ^ 0.000000000 

250.000000000 

260.000000000 

270.000000000 

280 . 000000000 

290 . 000000000 

300.000000000 
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INITIAL INPUT CONDITIONS 



ATMOSPHERIC INPUT PARAMETERS 

D (RELATIVE ATMOSPHERIC ROTATION) = 1.0000 

F107 (DECIMETRIC SOLAR FLUX) = 100.00 

FBAR (90 DAY AVG. F107) = 100.00 

AKP (GEOMAGNETIC INDEX) = 2.00Q. 

TVE (TIME OF VERNAL EQUINOX PASSAGE) = <:*3222 . 738200000 MODIFIED JULIAN DAYS 



SATELLITE PHYSICAL PARAMETERS 

CD (DRAG COEFFICIENT) = 2.20 UNITLESS 

A (CROSSECTIONAL AREA) = 0.000050000 SQUARE KILOMETERS 
SPI (MOTOR SPECIFIC IMPULSE) = 230.00 SECONDS 

W (INITIAL FUEL MASS) = 150.00 KILOGRAMS 

WT (INITIAL SATELLITE MASS) = 9100.00 KILOGRAMS 



NAVSPASUR ORBITAL ELEMENTS 

UJD (EPOCH) = ^4619.98716775 MEAN JULIAN DAYS 
ETU (MEAN ANOMALY) = 150.0000 DEGREES 

HO (MEAN RIGHT ASCENSION OF THE ASCENDING NODE) = 95.0000 DEGREES 

GO (MEAN ARGUMENT OF PERIGEE) = 200.0000 DEGREES 
XI (MEAN INCLINATION) = 95.0000 DEGREES 

RND (RATE OF CHANGE OF MEAN MOTION) = 0.000783912 REVOLUTTIONS PER DAY SQUARED 

ES (MEAN ECCENTRICITY) = 0.0100 UNITLESS 

ESD (ECCENTRICITY DECAY RATE) = -0.239590000 1/SECONDS 

AO (KAULA SEMI-MAJOR AXIS) = 1.06000 MEAN EARTH RADII 

ADOT (SEMI-MAJOR AXIS DECAY RATE) = -701.865200000 



CONTROL INPUT PARAMETERS 

MAXIMUM REVOLUTIONS PER ITERATION 5000 

NUMBER OF SCHEDULED ORBIT ADJUSTS TO SEMI-MAJOR AXIS 0 



L 



PROPELLANT LIFE ESTIMATOR 
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CDA= O.llOOOE-03 SQ KM 

SPECIFIC IMPULSE= 230.00 SEC 

INITIAL HEIGHT OF PROPELLANT= 150.00 KG 



F107= 100.00 

FBAR= 100.00 

TIME OF VERNAL EQUINOX= 0 . 43222738200000E+05 



MAXIMUM NUMBER OF ORBITAL REVOLUTIONS= 5000 



DATA PRINTED EVERY 100 ORBITAL REVOLUTIONS 



K02AI MEAN ELEMENT SET 

SEMIMAJOR AXIS= 0 . 675620541 077 90E+04 KM 
ECCENTRICITY= 0 . 965611329046 15E-02 
INCLINATION= 0 . 9499 996 48496 17 E+02 DEG 
ARGUMENT OF PERIGEE^ 0 . 19333742463643E+03 DEG 
ASCENDING NODE= 0 . 9500014966847 OE+02 DEG 
MEAN ANOMALY= 0 . 15661 152641309E+03 DEG 



REV NUMBER 


TIMECMJD) 


TIMECDAYS) 


PROPELLANT REMAININGCKG) 


100 


44626.32 


6.33 


145.28 


200 


44632.71 


12.72 


140.69 


300 


44639.10 


19.11 


136 .21 


400 


44645.49 


25.50 


131.83 


500 


44651.88 


31.90 


127.54 


600 


44658.27 


38.29 


123.37 


700 


44664.67 


44.68 


119.33 


800 


44671.06 


51.07 


115.45 


900 


44677.45 


57.46 


111.72 


1000 


44683.84 


63.86 


108.07 


1100 


44690.23 


70.25 


104.42 


1200 


44696.63 


76.64 


100.72 


1300 


44703.02 


33.03 


96.94 


1400 


44709.41 


89.42 


9.3.05 


1500 


44715.80 


95.82 


88.96 


1600 


44722.19 


102.21 


84.60 


1700 


44728.59 


108.60 


79.92 


1800 


44734.98 


114.99 


74.98 


1900 


44741.37 


121.33 


69.93 


2000 


44747.76 


127.77 


64.97 


2100 


44754.15 


134.17 


60.27 


2200 


44760.55 


140.56 


55.92 


2300 


44766.94 


146.95 


51.91 


2400 


44773.33 


153.34 


48.18 


2500 


44779.72 


159.73 


44.59 


2600 


44786.11 


166.13 


41.07 


2700 


44792.51 


172.52 


37.58 


2800 


44798.90 


178.91 


34.10 


2900 


44805.29 


185.30 


30.53 


3000 


44811.68 


191.69 


26 .78 


3100 


44818.07 


198.09 


22.74 


3200 


44324.47 


204.43 


18.37 


3300 


44830.86 


210.87 


13.66 


3400 


44337.25 


217.26 


8.68 


3500 


44843.64 


223.65 


3.58 


3600 


44850.03 


230.05 


-1.49 


3600 


230.05 




0.000000000 
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